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ABSTRACT 


A fully factorized two-dimensional Navier-Stokes flow solver has been 
developed and applied to the problem of predicting subsonic airfoil flutter in the 
light stall regime. The inviscid fluxes are evaluated with a central difference 
ADI scheme and fourth and second order numerical dissipation 1s used to obtain 
oscillation-free solutions. The performance of algebraic and one-equation 
turbulence models in predicting separated flow is explored for computing high 
Reynolds number steady flow and unsteady flows over an oscillating NACA 0012 
airfoil. Comparisons of the computed results with available experimental data 
indicate that even though the lift response is fairly well predicted, the 
computation of the pitching moment hysteresis loops is very sensitive to 
turbulence modeling. Results computed with several current models are in good 
agreement whenever the steady stall angle is exceeded only slightly. However, 
they fail to capture the vortex shedding process leading to the onset of stall 
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I. INTRODUCTION 


A. BACKGROUND 

The development of numerical solution methods for two-dimensional 
Navier-Stokes equations during the past few years provides new tools for the 
investigation and prediction of airfoil flows. Of great interest is the study of flow 
separation on airfoils in unsteady motion which is usually referred to as the 
dynamic stall phenomenon. McCroskey et al [Ref. 1| performed a series of 
careful experiments which serve as valuable benchmark data. More recently, 
Lorber and Carta [Ref. 2] contributed important additional experimental 
dynamic stall information for a Sikorsky airfoil. They also investigated incipient 
torsional stall flutter (Ref. 3] and found that small-amplitude airfoil oscillations 
near static stall may be unstable. 

For a pitching airfoil the instantaneous work done on the fluid by the airfoil 
due to its motion is the product of the pitching moment about the axis of rotation 
and the differential change in angle of attack. This product usually is a positive 
quantity. However, if the net work per cycle of oscillation were to become 
negative then the fluid would be doing work on the airfoil. Once the airfoil 
begins to extract energy from the freestream, the amplitude of the oscillation 
will grow and finally diverge. This condition is known as stall flutter. Normally, 
flutter of an airfoil is due to a combination of torsion and bending. However, in 
this case flutter is caused by a single degree of freedom oscillatory motion. 
During a cycle of oscillation the lift coefficient and the pitching moment 


coefficient plotted versus angle of attack produce a hysteresis loop. It is the 


pitching moment hvsteresis loop that provides an indication of incipient stall 
flutter. As shown by Carta and Niebanck [Ref. 4], clockwise pitching moment 
loops represent negative aerodynamic damping and therefore cause oscillations 
of a free airfoil to grow in amplitude, while counterclockwise loops cause such 
oscillations to decav. Hence torsional flutter will occur as soon as the area of 
the clockwise loop exceeds that of the counterclockwise loop. The 
aerodynamics of small-amplitude airfoil torsional oscillations near stall need to 
be investigated experimentally and computationally in order to examine the stall 


flutter mechanism. 


B. PURPOSE 

The first objective of this investigation is to test an unsteady, compressible 
Navier-Stokes code (Ns2.f) using an Alternating-Direction-Implicit (ADI) 
scheme based on the the Beam-Warming [Ref. 5] approximate factorization 
method to determine its ability to obtain realistic airfoil flow solutions for a 
variety of flow regimes is examined. The accuracy of the numerical solution is 
investigated by comparing the computed solutions with available experimental 
data. Test cases include steady-state flow solutions at various flow speeds and 
angles of attack, as well as, unsteady flow solutions over rapidly pitching and 
harmonically oscillating airfoils. In addition, the accuracy of the Navier-Stokes 
solutions are further explored by comparing each case with an unsteady, 
inviscid, incompressible, panel code (U2diif.f), and with a steady, 
incompressible, viscous/inviscid interaction code (Incompbl.f). The final and 
principal objective of this study is to determine the influence of mildly separated 
flow during part of a harmonic oscillation cycle on blade stability. Results are 


presented for NACA 0012 airfoils showing the influence of various parameters. 


Flow field details are also included in order to provide a better understanding of 


certain flow features that may lead to stall flutter. 


II. GOVERNING EQUATIONS 


The continuity equation, the momentum equation, and the energy equation 
must be solved simultaneously in order to obtain a flow solution of a 
compressible viscous fluid about a body. A complete derivation of these 
equation can be found in various texts, eg., Anderson [Ref. 6]. However, main 


steps of the derivations are presented in the following sections. 


A. CONTINUITY EQUATION 

The continuity equation is the result of applying the physical principle of the 
conservation of mass to a finite control volume fixed in space. Simply stated 
the net flow out of a control volume through its bounding surface must be equal 
to the time rate of change of the mass inside the control volume. For any 


arbitrary control volume, the continuity equation can be expressed as 








Ə p - 
——- V|pV |z0 
dt 2 | (1) 
which for a two-dimension Cartesian Coordinate system becomes 
op a 9 
(2) 


B. MOMENTUM EQUATIONS 
The equation for the conservation of momentum is obtained by applying 
Newton’s second law, which state that the net force acting on a fluid particle is 


equal to the time rate of change of linear momentum of the fluid particle. In 


Cartesian coordinates the momentum equations can be expressed as follows 


alpu) Aleu +p) 9(puw) 2r, дт, 
ex Ох z OX д (3) 
and in the y-direction. 
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Stress terms are as follows 
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C. ENERGY EQUATION 


The conservation law form of the energy equation is derived by applying 


the first law of thermodynamics (dE=dQ+dW) to a fluid particle. This leads to 


Æ д 

E ((e pa] [tE +p)w]= 
EE UE 1 о) 

Ox № (6) 


where the total energy per unit volume and the heat flux are given by 


Е = (е+ (27 јр, а=-К дТ/дх, (7) 


Неге К 15 the thermal conductivity. 


D. CONSERVATION FORM OF GOVERNING EQUATIONS 

The starting point for the numerical algorithm which is presented in the next 
section is the strong conservation law form of the two-dimensional Navier- 
Stokes equations. The non-dimensionalized vector form of the governing 


equations in conservation law form for a Cartesian coordinate system 1s: 
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Here the pressure is related to the conservation variables of q by the equation 


р=(у-1е- 0,5р(1? +?) (11) 


where y = 1.4 is the ratio of specific heats, and a is the local speed of sound 
given by, a2 2 y p / p. 

The density and the velocities are non-dimensionalized with the freestream 
density Pee and the free stream speed of sound ace, respectively. The total 
energy 1s normalized by р The time (t) scales as t* = t ao / c, where (c) 
is a characteristic length, such as the chord length. The Euler equations are 
obtained from Equation (8) by dropping the viscous terms. The strong form is 


chosen because it enables shock capturing. 


Ш. SOLUTION METHODS 


A. POTENTIAL FLOW METHOD (U2DIIF.F) 

For inviscid incompressible flow the conservation equations reduce to the 
Laplace equation which can be solved in a number of ways. We present here a 
brief outline of the widely used panel method for both steady and unsteady 
airfoil flow. 

I: Panel Method for Steady Airfoil Flows 

The frame of reference for the formulation of the steady flow problem 
is a fixed (x,y) coordinate system located at the leading edge of the airfoil. The 
freestream velocity is represented by +Voo and the x-axis passes from the 
leading edge through the trailing edge of the airfoil. 

In order to formulate a method for computing the flow around an 
airfoil, the airfoil surface is divided into (n) straight line segments or panels. The 
end points of the panels are called nodes, and since there are (n) panels there 
are then (n+1) nodes. Complex airfoil geometries can therefore be modeled 
with a greater number of nodes and panels. The trailing edge panel is the first 
panel, and the next panel continues on the lower surface in a clockwise manner 
until the ni^ panel is reached at the trailing edge of the upper surface. 

Now that the frame of reference has been specified, two additional 
vectors must be defined, the unit normal vector nj, which is always 
perpendicular to the (ith ) panel and directed outward from the airfoil surface 
and the unit tangential vector tj, which is parallel to the (it!) panel and is 


directed from the (n) node to the (n1) node. 


Two types of singularities are sufficient to model lifting airfoil flows. 
U2diif.f places uniform source distributions qj and a uniform vorticity distribution 
y on each of the j-panels. Both singularities satisfy Laplace's equation, a linear 
homogeneous second order partial differential equation. Since the solutions to 
linear PDE’s can be superimposed on each other, a simple flow can be added to 
another simple flow and so on, until a very complicated flow field is created. 

The flow field around an airfoil, represented by the velocity potential, 
can be constructed from the potential of the freestream flow added to the 


velocity potential of source and vorticity distributions, hence 


P = @. + @, + 
Фф. + Ф, T 9, (12) 
where, 
Q, = У. +(Х соза ғу sin G) 
(13) 
(p, = | WS, rds 
(15) 
| Хш 
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(15) 
summation of the three potentials yields 
P = V (x cosa y sin æ)+ Í W rds- Í 76) 0 ds 
2л 2л (16) 


Equation (16) can now be evaluated at any point in the flow field by evaluating 


the integrals along the airfoil contour (s), where the flow field point is located at 


a distance (r) at an angle (0) measured from a point on the airfoil. Introducing 
n panels and then summing the contributions of each panel yields 
qd; 20 


In Le jas 
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(17) 
Once the boundary conditions are defined the system of (n) equations can be 
solved for (n+1) unknowns. 

Next, 1t 1s useful to introduce the concept of influence coefficients. An 
influence бее is the velocity induced at a point in the flow field (field 
point) by a unit strength singularity, source or vorticy, placed anywhere within 
this field. U2diif.f places these points on each panel. A detailed description of 
the use of influence coefficients is found in Teng [Ref. 7]. 

For the steady flow problem the influence coefficients are given in 
Table 1. 

TABLE 1. INFLUENCE COEFFICIENTS 





Normal component induced at ith control point by 
unit source distribution on j^ panel. 
Tangential component induced at it control point 


by unit source distribution on jt” panel. 


Normal component induced at ith control point by 


unit vorticity distribution on jt? panel. 


Tangential component induced at ith control point 





by unit vorticity distribution on jth panel. 
Angle of control point (i) and panel (j) measured 


from x-axis. 
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There are only two boundary conditions that must be satisfied in order 
to solve the lifting airfoil flow problem. The first boundary condition is the 
condition that the flow remains tangent to the airfoil surface and the second 
condition is that the pressures on the upper and lower trailing edge panels must 
be equal. This condition is known as the Kutta condition. Pressure is related to 
velocity through Bernoulli's equation for steady potential flow, which allows the 
Kutta condition to be expressed as 

(venie) да „дап ын | 
ower upper (18) 
The flow tangency condition is expressed as 


(vm) -0 1=12,.., п 


1 


(19) 
Equation (18), the Kutta condition, using the influence coefficient concept is 


expressed as follows 


M^] -Y2iBi - V.cos(a -6,) - 
йн = 


2 (Аға | уу [в: |+ V. cos(a - Ө.) 
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Equation (19) the flow tangency condition, becomes, using the influence 


coefficient form 


У [Ага |+ уу [в |+ V. sin(a -0) и 
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These equations can then be written in. the following matrix form 


ЕТ SI Па, | | bi | 
| DLE 
س‎ 
: q, b 
me ко Hn г м (22) 


and solved with the method of Gauss Elimination with Partial Pivoting. 
2. Unsteady Numerical Formulation 

A brief and concise description of the unsteady numerical formulation 
is found in Krainer [Ref. 8]. To model unsteady flow around an airfoil using the 
panel method, N unknown source strengths and one unknown vorticity strength 
are located along the airfoil, a wake panel of unknown vorticity strength, length 
and orientation is attached to the trailing edge of the airfoil. This makes a total 
of N+4 unknowns. The flow tangency requirement at each of the N panels 
represents a system of N equations. Using influence coefficients and using the 
subscript k to count time, the tangential flow condition of the itf element is 


written as follows 


n 


Ума) Би, ув + 12 -(0(01-У(0)-000(х1- ЛЕ 


j=l 


к-1 
зо (Bi), 5 28 55 ээ Е" Г) | mU це n 
m=] 


k 


(23) 
where Voo 1S the mean velocity, (U(t)i + V(t)j) is a time dependent translational 
velocity and Q(t) is a rotational velocity. The vectors i and j are in the airfoil 


fixed coordinate system. 
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The Kutta condition is a single equation that determines the circulation 
around the airfoil. In this case the pressures at the upper and lower surface are 
equated at the midpoints of panels adjacent to the trailing edge (1.e., the first and 


last panel) 


The Helmholtz theorem is another single equation that provides a 
relation for the strength of the vorticity distribution along the wake element. 
The Helmholtz theorem states that any change in circulation around an airfoil 
must be countered by a change in vorticity in the wake of equal magnitude but 


of opposite sign. This theorem can be written as follows 


ПИ EIE) 
(252 


where /is the length of each individual wake element and yits vorticity 
strength. Therefore, the vorticity in the wake element is equal to the negative 
change in circulation around the airfoil with respect to the k-1 timestep. 

To acquire the two additional equations required to solve for the N+4 
unknowns, assumptions about the geometry of the wake panel are made. First, 
the wake panel is oriented in the direction of the local resultant velocity at its 


midpoint, as viewed in the (moving) airfoil frame of reference. 


(v), 


tanO, = 
2 (26) 
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(VYw)k and (vXy)& are the x- and y-velocity components at the midpoint of 
each wake panel. Second, the length of the wake panel is assumed 
proportional to the magnitude of the local resultant velocity at its midpoint and 


to the timestep. 


шэн (4), | + (уз М te -t i] 


The nonlinearities in the Kutta condition and in the wake panel 


(27) 


assumptions necessitate an iterative solution procedure. 


B. VISCOUS/INVISCID INTERACTION METHOD 
(INCOMPBL.F) 


The theory and numerical methods presented in this section are taken from 
the work of Cebeci and Bradshaw [Ref. 9] and the investigations performed by 
Krainer [Ref. 8] and Snir [Ref. 10] AS stated in the introduction the results 
from the Navier-Stokes code (Ns2.f) were compared with a viscous/inviscid 
interaction method code (Incompbl.f) made available by Cebeci. An 
abbreviated description of this method is given in this section. For a complete 
description see the original publications, References 8 and 9. 

The governing principle behind the viscous/inviscid interaction method 15 an 
approximation to the Navier-Stokes equation that allows a flowfield to be 
divided into an inner viscous region and an outer inviscid region when the 
Reynolds number is sufficiently large. This concept of a thin viscous layer near 
the surface of a body and an outer region where these viscous effects are small 
compared to the inertial effects is known as the boundary layer theory. Unlike 


the Laplace equation, which is the governing equation for potential flow, the 
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thin shear layer equations are nonlinear 





СОН ду =0 
dx ду 
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(25) 
with the boundary conditions 

Zo mo 
у=оо, и-це(х) (29) 


The direct boundary layer method solution to the boundary layer equations 
by the Keller Box Method involves four steps. First the boundary layer 
equations are transformed into a system of first order differential equations. 
Next the Keller box method is used to approximate the first order differential 
equations by simple centered differences and two-point averages, using values 
at the corners of one difference molecule only. Newton's method is used to 
linearize the resulting algebraic equation. The Keller block elimination method 
is then used to solve the resulting block tridiagonal system. The procedure 
described above does not provide solutions to flows that are separated or have 
regions of reverse flow. 

1. Interaction Method 

The interactive boundary layer method provides a special coupling 
between the inner viscous flow and the outer inviscid flow, which enables 


reverse and separated flows to be calculated. In such areas, the external 
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velocity is substantially changed by the viscous effects and can no longer be 
considered as a known boundary condition for the boundary layer flow. 

The general approach to the solution is the same as for the direct 
method with modifications. Since the outer flow is unknown, the velocity at the 
edge of the boundary layer is given by the interaction law and is written as 


11 484 ғу 4 
CO а) 5: 





(30) 
where ue(x,ye) is the total velocity at the edge of the boundary layer, uej(x) 15 
the velocity computed by the inviscid panel method, 6* is the displacement 
thickness, and the integral term is known as the Hilbert integral. 

The following transformations are used to transform the boundary 


layer equation into a system of first order differential equations 
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(31) 
The boundary layer equation takes the form 
PSU RS N 
(БУ J+ xw = 12:523 
У (32) 
with boundary conditions: 
n=0, U(x,0)=0, f(x,0)=0, 
N=Ne, U(x,ne)=W(x,ne) (33) 
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The velocity at the edge of the boundary layer now becomes 


we) = ی‎ а lw ern. Ел 25 


The finite difference box method is used to solve the equations, in the 


(38) 


same way it was used for the direct case, but with two additions. First in areas 
of flow reversal the udu/dx is omitted to assure stable integration. And second, 
the edge velocity 1s approximated by the relation presented in the next section. 

By using central differencing to approximate the differential equations, 
a system of nonlinear algebraic equations is obtained for the unknown variables 
(which are f, U, V, and W). To solve the system of equations, the system is 
linearized by the Newton iterative procedure, and the resulting linear system 18 
solved for the new unknown variables which are the increments of, óU, ОУ, апа 
OW. 

The solution procedure is repeated until the change in the increments 
is negligible compared to the preceding iteration. The iterative process is 
performed again at the next downstream station. 

2. Interactive Model 

The interactive model is used to couple the boundary layer to the 
external flow. It is needed in areas where strong interaction occurs, and both 
the boundary layer and the outer flow must be solved simultaneously. The 
external velocity is assumed to consist of a potential flow term and a correction 


term due to viscous effects. [Ue(x )=Ue](x )+Ue§(X)] 
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The normal velocities at the surface of the airfoil, induced by sources 
distributed on the surface, displace the streamlines from the surface in the same 
way that the actual boundary layer displaces them. 

Several assumptions are made in order to express the correction term 
in the form of the Hilbert integral. First, the surface is approximated to be a flat 
plate, and the normal velocity is then half the local source strength G(x). 
Second, the inviscid velocity does not change across the boundary layer. 
Therefore, the local horizontal velocity induced by the source distribution, is the 


correction term to the inviscid velocity, and can be represented by the Hilbert 





integral 
š 1а раа 
пе = ЛС” E (35) 


The integration is carried out on all the sources on the surface. The Hilbert 


integral is then approximated by a finite series 


1 d 
и,6 )- cafu.) 





where Cik 1S a matrix of interaction coefficients which are functions of the 
geometry only. 

Since the computation of ue§ involves values of 6* downstream of the 
current x location, which are not known yet, these terms are taken from the 
previous iteration using a relaxation formula. 

3. Turbulence Model 

The turbulence model used in the code is the Cebeci-Smith model 

[Ref. 11] and Michel’s method is used to predict the transition from laminar to 


turbulent flow. 
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C. EULER AND NAVIER-STOKES CODE (NS2.F) 
In order to facilitate the numerical implementation for arbitrary complex 
flow domains, the Navier-Stokes equations are transformed to a generalized 


coordinate system, (6,6), using the transformations 


8 = 6л) 

44452 (37) 
The grid spacing in the curvilinear space is uniform and of unit length so that 
unweighted differencing schemes can be employed for the numerical 
implementation. The grid points in the Cartesian system, referred to as the 
physical domain, correspond via a one to one relationship to the points of the 
curvilinear transformed system, referred to as the computational domain. 
Singularities of the transformation may occur on the computational boundaries. 
Such singularities occur for grids with multiple connected regions. The 
transformation of the governing equations from the physical domain to the 
computational domain is obtained in most cases by numerically evaluating the 
metrics and the Jacobian of the transformation. The derivatives with respect to 
the x, y variables can be expressed in terms of the new variables by the chain 


rule. 


Жэ 1. 
E 
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The metrics and the Jacobian of the transformation are given by 


cx J 2 Cx = -J zg, (39) 
62 = -] T Cx =J xg, (40) 
| = = 

(2; г (41) 


The Navier-Stokes equations written in generalized curvilinear coordinates 
retain the strong conservation law form expressed by Equation (8). The 


governing equation for generalized curvilinear coordinates is 


^ ӘС 
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The conservation variable vector d and the inviscid fluxes Е апа С in the 





~ |" 


transformed space are given by 


Гр 1 [ pU 1 [oW | 
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п. Зи m am ри Єр | ~ J! pwW+¢,p 
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where U and W are the contravariant velocity components along (ће < апа С 


directions respectively, given by: 


U=6, tutiw, W= tut w (44) 


The viscous flux terms are transformed as 


^ 
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with the stress terms from Equation (10) expressed in terms of the transformed 


variables 5, С аз 
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(46) 
When the thin layer approximation [Ref. 12] to the two-dimensional 
conservation law form of the governing equations is applied, they take the 


following form 





24 2Е 96 _ 1 [06, | 
dt & de Rel a | (47) 
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Гр | | pU | | pW | 
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The viscous flux term is transformed as 


( 0 \ 

ё Е. ити, + (u / 3)m, Ç, | 
" J| H m,w,+ * (u/3)m, б, 

E m,m, + (u/3)m, m, | 


(49) 
Here, 
| НС. 
i, = Ue +Ç,W, 
1 2 | 2a2! 
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т, = буй+ бу sm 


1. Numerical Grids 

In order to compute flow solutions of partial differential equations 
(such as the Navier-Stokes Equations) with finite differences, a discretized 
version of the physical domain must be generated. During the course of this 
investigation, several different NACA 0012 airfoil C-type grids were generated 
using two different grid generation methods. A grid generation code call 
GRAPE that uses the Poisson differential equation was used in the early stages 
of investigation. The inviscid and viscous grids shown in Figures 1 and 2 were 
generated from an algebraic grid generation code. Note the finer resolution of 
the viscous grid near the airfoil surface. This resolution is required to 
accurately represent the boundary layer. The flow solutions presented in 


Chapter IV were computed by using the grids shown in Figures 1 and 2. The 
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advantage to using an algebraic grid , vice a PDE method or a conformal 
mapping method, is that it is computationally efficient. 

A simple description of the algebraic grid generation method is offered. 
First the airfoil surface coordinates are unwrapped to form a simple curve, by 
separating the lower trailing edge from the upper trailing edge. This curve 15 
the starting point for the generation of the computational domain. Lines are 
drawn normal to the curve and grid points are then located at intervals on these 
lines. Lines may be clustered near the leading or trailing edge of the airfoil; or 
at some other area of interest, eg., a shock location. The resolution of the 
points normal to the surface is achieved by the use of an algebraic function that 
provides uniform stretching normal to the body surface. Once the grid is 
created in the computation domain it is then transformed back into the physical 


domain using inverse transformations. 
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Figure l. NACA 0012 Inviscid 151x64 C-Grid 
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2. The Numerical Scheme 
The integration in time is performed implicitly with a second order 


accurate Scheme of the form 


n nel n_Ato/, n)\, Ata з | 
Aq =q q = == =- (ла += E ноја" гаа 44 


Using this trapezoidal rule differencing scheme to approximate the time 
derivative of the conservative dependent variable q vector, the unknown 


conservative variables at the n+1 timestep are given by 


му на Ea (а) (52) 


The spatial derivatives of the governing equations are discretized using central 
differences, and the right hand side of Equation (52) becomes 
n At ^ ^ Anl 
Aq --A Af «5,6 G) +(5,F +6,G) | 


2 (33) 


The nonlinear terms F and G at the n+1 timestep are linearized as follows 


A 


раз! P n+1 [A2)= F" + A" Aq"! +0 (42) 





(54) 
where A” is the flux Jacobian matrix given in Appendix A. Substitution of the 


linearized form of the flux vectors in Equation (54) yields 


lij - (8,1, б, мед - Afi - 5.6", } (55) 
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The two-dimensional operator on the left hand side of Equation (55) is 


approximately factorized as follows: 


11 865, ш 205,87, (аа = (fn, 3,67] 
(56) 


On the right hand side of Equation (56) an explicit dissipation term Dexpl, is 
added for numerical stability. In order to obtain an oscillation free solution a 
second order. dissipation term Dimpl 15 added to the left hand side spatial 
operators. Thus, the complete discretized form of the governing equations 


becomes 


A A + 
Ше (8,A, + A іш Ер) ва = 


(57) 
I-A [&F^, = 5,6 л m оо = (RHS) (58) 
This equation is solved by performing two sweeps as follows: 
ІШ» 6 А“ + (Рај Әм ix =(RHS)" 
2 ¢ Ik шар! (59) 
„n+l 
fu t (6,6 B D...) J = Ла ак (60) 


During each sweep the following linear block tridiagonal system of equations 15 


solved 


а, АФ +6, АФ + с. Aq... = f, (61) 


2] 


where Ag = Aq +! and f = (RHS)" when the Š direction sweep is performed, 
and Aq 2 Aqn*l and f - (RHS)*" when the € direction sweep is performed. 


The block matrices have the form 





206 Ji 
а; к DE: . E (62) 
i+] ,k 
b, - (1-26, A Jul 
(63) 
At F. Ше 
Cik = Aix б арія ] (64) 
i-lk 


3. Boundary Conditions 
All flows were computed at subsonic free stream speeds. For 
subsonic inflow and outflow boundaries the flow variables are evaluated using 
zero order Riemann invariant extrapolation. At the inflow boundary there is one 
incoming and three outgoing characteristics, therefore, three variables, density 
(p), normal velocity (wi), and pressure (p) are specified and the fourth variable 
axial velocity (ui), is extrapolated from the interior. The inflow boundary 


conditions are given by 


a iD (y -1) 
pi 3 , 5| 121221 d, = -R;) 


75 4 
и; = (Rt +В; )/2 
W, = w. 


2 
Di = (2) 
4 


(65) 
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where R*1, R72 are the incoming and outgoing Riemann invariants given by 
R;-u,-42a,/(y-1, Rj =u,—-2a,/(y-1) 

(66) 
At the outflow boundary there is one incoming and three outgoing 
characteristics and only one quantity, pressure, is specified while the others are 
extrapolated from the interior. For the density and the normal velocity, simple 
first-order extrapolation is used, and the axial outflow velocity is obtained from 
the zero order outgoing Riemann invariant. The outflow boundary conditions 


are given by 


А = р› 
MR 2D di = МУР. Ир, 
\ = \> 
P, Р, 


(67) 
On the body surface the nonslip condition is applied for the velocities. The 
density and pressure are obtained from the interior by extrapolation. For the C- 
type grids used in this study averaging of the flow variable at the wake cut is 
used. 
4. Turbulence Models 
a. Baldwin-Lomax Turbulence Model 

The Baldwin-Lomax model [Ref. 12] is a two -layer, inner and 
outer eddy viscosity model for the computation of two- and three-dimensional 
flows. Patterned after the Cebeci-Smith model [Ref. 11], it incorporates a 
modification that bypasses the need for finding the edge of the boundary layer. 


This is achieved by introducing the vorticity in place of the boundary layer 
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thickness. The turbulent effect is simulated by an inner eddy viscosity, given 


by: 


2 
| == < 
(LL, ree pl k»| y Sue (68) 
where ycrossover is the smallest value of y at which the values of the inner and 
Outer eddy viscosities are the same and the Prandtl mixing length 115 
яти о 
The Ж ) 
(69) 


The magnitude of the vorticity in two-dimensions lol and y* are defined as 





follows 


"РЕ 21) yt = Petey LN Pets Y 
ду ox } and Hy LL, (70) 


where A* is a constant. The subscript w denotes values at the wall or airfoil 
surface in this case. 

The outer eddy viscosity is given by the following expression 

(H outer = KG E в ES y (] l) 
where к апа Ccp are constants and F(y)wake=Ymax Fmax for boundary layers апа 
FO) wake = wee (UpIF2/Fmax ) for wakes and separated boundary layers. 


The quantities Ymax and Fmax are determined from the expression 


a "a 
Б(у) =у [о ! (3 1 ^ 





30 


The exponential term of Equation (72) is set to zero for wakes. Fmax is the 
maximum value of F(y) that occurs in a profile and y,,,, is the value of y at 


which it occurs. F(y)kygng is the Klebanoff intermittency factor given by 


| 3 
Fa над E | "i 


The quantity Upp is the difference between Umax=Uy=ymax and the minimum 


total velocity in the profile 


Upr 7 (Vu? v? ] (К+) 


Y= You У У ш (74) 
The second term in Upir is taken to by zero except in the case of wakes. 

The constants were determined to achieve agreement with 
Reference 11 and can be found in the original document. This model 
completely models turbulent flow over an airfoil, it is relatively easy to code, it 
does not degrade the solution convergence, and it is computationally efficient. 

b. Johnson-King Turbulence Model 
The Johnson-King model [Ref. 13] takes into account the 


convective and diffusive effects on the Reynolds shear stress -u'’w’ in the 


streamwise direction. The eddy viscosity is given by 


| С | 
V = V. —ехр — 
| i Vio | (75) 


where vt , Vto describe the eddy viscosity variation in the inner and outer part of 


the boundary layer. The inner eddy viscosity is computed as 
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ә 2 TT eT айн + 
и, = D Ky (cu Vies where 0? =1-е А | where the constant A*-15. The 
outer viscosity is given by 
Vto*o(x)[0. 1 68U46*y] (76) 


where y is the Klebanoff intermittency function y-[145.5(y/8)9]-! and o(x) is the 

solution of the ordinary differential equation which describes the development of 
-u'w'lmax along the path of maximum shear stress. This model accounts for the 
effects of convection and diffusion on the Reynolds stress development through 


the solution of the following ODE 











| ] 
Z ай 1-2) Cat La - Ч, | 
Ox 20.) 29 17-2) У, ед | 


| 9 


here Cqif and aj are modeling constants, um is the maximum average mean 


(77) 
velocity and g=[-u’w’lmax]-!/2, geq7[-u' w'Imax,eg]- 1? where Lg is the 
dissipation length evaluated as 


Г. =040 у, /5<0225 


І, «009 620225 
бі у Yd] (78) 
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The equilibrium shear stress geg in Equation (77) is determined from the 


following equilibrium eddy viscosity distribution 





[ E ? ) 
M I скр = 
eq "| oes | 


У ед = Реку (CU... 


Vi oq 70(x)0168U,6 y | 


19.64 (79) 
An implicit Euler method is used for the numerical solution of Equation (77), 


and the maximum shear stress at each iteration level is updated as follows 


n+l 


ск)" зах) t2- 
К (80) 
Solutions with the Johnson-King turbulence model were obtained as follows. 
First a convergent solution using the Baldwin-Lomax turbulence model for the 
entire flow field was computed. Then the Johnson-King model was applied only 
to the upper part of the airfoil. To initiate the solution o(x) in Equation (76) is 
set equal to unity and it is allowed to change according to Equation (80) until 
the final solution is obtained. 
c. Algebraic RNG-based Turbulence Model 

Recently an algebraic eddy viscosity, as well as a two-equation K- 
€ model based on Renormalization Group (RNG) Theory of turbulence [Ref. 14] 
were proposed for the closure of the Reynolds-averaged Navier-Stokes 
equations. The algebraic model, although free from uncertainties related to the 
experimental determination of empirical modeling constants, still requires 


specifications at an integral length-scale of turbulence which reduces the 
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generalities of the model. Here the integral scale is assumed proportional to the 
distance from the wall and the eddy viscosity for the RNG-based algebraic 
turbulence model is obtained as in Martinelli [Ref. 15] from the following 


formula 
г | Т 
^ d 3 
а |l 1 \ 
нэрэ 24388 "€ 
(81) 
where v=vt+v], H is the Heavyside step function and @ is the dissipation 
function $—tij ( dui / 0x; ). The RNG-based turbulence model is applied only for 


the suction surface separated flow region while the pressure side and the wake 


regions are computed with the Baldwin-Lomax model. 
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IV. RESULTS AND DISCUSSION 


A. METHOD OF INVESTIGATION 

As stated in the introduction this investigation is divided into three parts; 
validation of the Navier-Stokes code (Ns2.f) through comparison with 
experimental data and a well tested inviscid unsteady Panel method code as 
well as a steady viscous/inviscid code. Finally evaluation of the effects of 
reduced frequency, mean angle, amplitude, Mach number, and Reynolds 
number on unsteady flow behavior are examined. Also, during the course of 
this study it became necessary to investigate the effects of several turbulence 


models on the results computed for the harmonically oscillating cases. 
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The following table lists the cases studied during the course of this 
investigation: 


TABLE 2. TEST CASES 


Type of Motion 
1 


Steady State M=0.7, a=1.4°, Re=9x 106 
(Transonic) 


Steady State 













М=0.799, а=2.26“, Ее=9х 106 





(Transonic) 


Steady State 
Steady State 


Rapidly Pitching Airfoil 





М=0.3, Ее=3х100 






М=0.3. Ве=4х106 







M=0.3, k=0.01272, Re=2.7x106 









(Ramp) 






Harmonically Oscillating QG(t)2AQ XA ] 'sin(OX) 







(Multiple Parameters) 






o:(t)=13°2.5°sin(wt), k=0.2, Re24x106 





Harmonically Oscillating 








(Small Amplitude) 





o:(t)=13°42.5°sin(wt), k=0.1, Re=2x100 





Harmonically Oscillating 


(Small Amplitude) 


LE Harmonically Oscillating 


Harmonically Oscillating С(0-1055.5 (00), К-0.1, Ёе-4х100 
Harmonically Oscillating o(t)=11°$5°sin(at), k=0.1, Re=4x 10° 






o.(t)=9°+5°sin(at), k=0.2, Re=4x 10 












1. Steady-State 

Flow solutions can be obtained either by rotating the grid or by rotating 
the flow. Cases 1 and 2 solutions were obtained by rotating the grid. 

First the procedure was used for the computation of steady flow over 
the NACA 0012 airfoil. Steady-state solutions were obtained for subsonic and 
transonic flow regime. Flow solutions for Case 1 and 2 (Table 2) were 
computed by rotating the direction of the flow relative to the stationary grid. 
However, in Cases 3 and 4 the grid was rotated to the desired angle using a 
simple Fortran code called rotgr.f and the free stream flow remained parallel to 
the x-axis. The program also translates the grid to a desired pivot point. In this 
thesis the pivot point was always chosen to be the quarter chord point. Either 
method produced the same solution. The steady-state solutions tested the 
accuracy of Ns2.f and provided initial solutions for the unsteady cases. Ns2.f 
was also compared to an existing steady, incompressible, viscous/inviscid 
method (Incompbl.f) obtained from Dr. T. Cebeci, California State University at 
Long Beach. 

a. Case 1. M=0.7, a=1.4°, Re=9x106 

The flow conditions for the first test case are M.=0.7, at an angle 
of attack of 1.4 degrees, and a Reynolds number of nine million. Explicit 
boundary conditions were used and the Baldwin-Lomax turbulence model was 
selected. The solution was obtained on a 161x64 point grid. This case was 
chosen because the solution converged rapidly. Figure 3 shows good 
agreement between the inviscid (Euler) pressure distribution and the measured 
values of Harris [Ref. 16]. Next, a Navier-Stokes solution was obtained whose 


convergence history is given in Figure 4. The solution was carried out to 2500 


ow 


iteration at a timestep of approximately 0.08. The pressure distribution 
computed by Ns2.f (Figure 5) is also in good agreement with the measured 
values. However, the suction peak is more accurately predicted by the inviscid 
solution. For a free stream Mach number (Mo) of 0.7 a maximum Mach 
number of 1.1 (Figure 6) is obtained on the upper surface near the leading edge. 
The density contour plot, Figure 7, shows the smooth contours indicative of a 
fully converged solution. The velocity field, Figure 8, shows the boundary layer 
thickness, represented by the thick dark line. The sonic line is represented by 


the line extending from 7% chord to 20% chord. 
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Figure 3. Pressure Distribution (Inviscid), M=0.7, a=1.4°, Ке=9х106 
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Figure 5. Pressure Distribution , M=0.7, a=1.4°, Re=9x106 


41 


MACH NUMBER 


0.700 MACH 
1.40 DEG ALPHA 
CONTOUR LEVELS ыш. 
Шор ТІМЕ 
0.00000 
0.05000 16 кб4 GRID 
0.10000 
0.15808 
0 806; (c) 
0.55000 
0.90000 
95060 s 
; 00000 К 
1.05000 ORT ыы 
0000 "ИК cede F 
Figure 6. Mach Contour, M=0.7, a=1.4°, Re=3x106 


IP ESSET 


0.700 MACH 
140DEG АІРНА 
9.00x10**6 Re 
| 16. TIME 

0.62000 TIME 

0.64000 

0.66000 

0.68000 

944000 

2000 


гээх 
ЫГ 2 


` Л. #2 
. "P 7 
%%442 “7, 


š. $ “we . 4... 
.. Ксы” , 


= ر‎ 
y Ы SX 
ees EA) 
з? ROE 


0 2" 
"ode Bo 


ОКО ЕЕ я 
10000 
12000 
14000 
16000 
16000 
20000 
22000 
2201) 
1.26000 


si 


paed «а-ға ь---» Ь«--а калай paad tinand лола 


Figure 7. Density Contour, M=0.7, a=1.4°, Re-9x109 
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Figure 8. Velocity Field, M=0.7, a=1.4°, Re=9x106 


b. Case 2. М=0.7 99, а=2.26°, Ке=9х106 

The next steady-state case investigated was at a Mach number 
equal to 0.799, at an angle of attack of 2.26 degrees and a Reynolds number of 
nine million. As in Case 1, an Euler solution was first obtained and compared to 
measured data as well as to a viscous solution obtained with an upwind scheme 
using the Johnson-King turbulence model. As seen in Figure 9, the pressure 
distribution for the inviscid solution is in fair agreement with the measured data 
until the shock is reached at the 50% chord point. After this point the computed 
solution fails to predict the pressure jump due to the presence of a shock at 
mid-chord. Note that the upwind scheme in combination with the Johnson-King 
turbulence model accurately predicts the location of the shock, even though the 
magnitude of the pressure gradient is off slightly. Next, a viscous solution was 
computed by Ns2.f with the Baldwin-Lomax turbulence model. Fourth order 
dissipation was added to produce a smooth flow solution across the weak shock 
located at mid chord. The solution was continued for 6000 iterations. The 
spike in the convergence history (Figure 10) was due to the addition of more 
fourth order dissipation. It is seen that the shock location was not predicted. 
At this angle of attack and Mach number the flow starts to separate and the 
Baldwin-Lomax turbulence model fails to model this separated flow region. 
Figure 11 shows the viscous pressure distribution comparison. The sonic line is 
shown on the Mach contour plot of Figure 12. A maximum Mach Number of 
1.45 is reached in this supersonic flow region. Note the steep Mach gradient 
(Figure 12) at approximately 6096 chord. Figures 13 and 14 show the pressure 
contour and density contour, respectively. Again the steep density gradient 


indicates the shock location at approximately 60% chord vice the measured 
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location at 5096 chord. Finally, the velocity field and sonic line are shown in 


Figure 15. 
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Figure 15. Velocity Field, M=0.799, и=2.26°, Ве=9х106 
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c. Case 3. М-0.3, Ке-3х106 

In order to produce C] vs œ and Cm vys @ curves, solutions were 
computed. faga? a-4 9-607025 7029 7 o0 50 IMEEM 
a=14", and a=15°. However, only results for a=4°, a=11°, and a=14", are 
presented. 

Figure 16 shows the convergence history of Ns2.f for a 161x64 
algebraic C-type grid. The Courant-Friedrichs-Levy (CFL) condition was 
approximately 2000-2500 which corresponded to a timestep between 0.005 and 
0.01. A comparison of the computed pressure distributions at a=4° (Figure 17) 
shows good agreement between the Navier-Stokes calculations using the 
Baldwin-Lomax turbulence model and the viscous/inviscid code predictions. 
The measured suction peak is somewhat higher than predicted by the codes. 
Figures 18 and 19 show fair skin friction agreement. Note that the Navier- 
Stokes code assumes turbulent flow over the entire airfoil, while the 
viscous/inviscid interaction code uses Michel's method to predict the 
laminar/turbulent transition. As seen in Figure 18, the transition point for the 
upper surface is at 10% chord. A Mach contour plot is shown in Figure 20 and 
a density contour plot with the normalized stagnation pressure (0.98) contour 
overlayed in Figure 21. A maximum Mach number of 0.48 was attained. The 
density contours are smooth with steep gradients at the leading edge. 

In Figures 22 through 29 similar results are presented for angles of 
attack of 11 and 14 degrees. It is seen that the agreement between the two 
codes and the measured pressure distributions is quite satisfactory, although the 


suction peak tends to be underpredicted by both codes. 
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Figures 30 and 31 show the Cj vs @ and Cm vs @ curves, 
respectively. Stall is seen to occur at approximately @=13.5°. The Navier- 
Stokes code reproduces the experimental lift and moment values up to the 
maximum lift value quite well, whereas the viscous/inviscid code shows greater 


deviation for angles of attack exceeding about eight degrees. 
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d. Case 4. M=0.3, Re=4x106 
Figure 32 shows a comparison of the computed lift curve with the 
experimental data at a Reynolds number of four million. It can be seen that the 
computed maximum lift angle and hence the static stall angle is 13.5° degrees. 
This agrees quite well with the experimentally measured stall angle of 
McCroskey et al [Ref. 17]. The computed lift and pitching moment are in fair 


agreement for small angles of attack, but start to deviate at larger angles. 
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M = 0.3, Re = 4 x 108 Baldwin-Lomax Turbulence Model 


Experiment: McCroskey 





2.0 — 

1.8 — 161 « 64 C-Grid 

- | « Measured 

| —— Computed 
aj 
Cos ds 
S 1.2 — 
9 | 
= 1.0 — 
о | 
52 
= | 
= .6 E 

4 = 22 

| 
21-25 
[7 | | | | | 
0 2 4 6 8 10 12 14 16 18 20 
Angle of attack, deg 

.20 = 

" | 161 « 64 C-Grid 

| | e Measured 
з 10 - "ЗУ Computed 
TE 
= Da 
z е ее 
£ O-—-*Oeacmassctettih cs 
zoo é п ыт 
Е-.05 — рс 
2 | Фо ое 
= 


4 
А 
e 


| 1 
N — 

e сл 

--т- 


0 2 4 6 8 10 12 14 16 18 20 
Angle of attack, deg 


Figure 32. C| vs. @ and Cm vs. a, M=0.3, Re=4x106 


73 


2. Rapidly Pitching Airfoil 

Next, AGARD Case 6 from Landon [Ref. 18] was investigated using 
Ns2.f and U2duf.f. The computed results are compared to the measured data. 
The flow condition for this case corresponds to a freestream Mach number 
equal to 0.3 at a non-dimensional pitch rate equal to 0.01272 and a Reynolds 
number of 2.7 million. Flow solutions for a rapidly pitching airfoil were 
computed by pitching the airfoil about the quarter chord at a constant rate from 
zero degrees angle of attack to a final angle of attack of 15.54°. Figure 33 


shows a sketch of the motion produced by a rapidly pitching airfoil. 


15.547 


Angle of Attack 





-3500 0 4000 


Iterations 


Figure 33. Rapidly Pitching Airfoil Motion | 


First a steady-state solution at zero degrees angle of attack was well 
converged to 3800 iterations. Then the iteration counter was reset to zero and 
the airfoil was rapidly pitched up to the final angle. The reduced frequency k is 
given by kza c/ 2 Ue,, where a is the pitch rate. In terms of non-dimensional 


quantities &(t)-2M.k and a(t)= a9+a1(t/ T), where (0 (Ao), is the initial angle 
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of attack «=0°, «1 (A1), is the final angle of attack &=15.54°, and Т is the time 
required for the airfoil to complete the pitching motion. 
a. Case 5. Mz0.3, kz0.01272, Ке-2.7х106 

Figure 34 shows the comparison of the computed pressure 
distribution for the ramp-type unsteady motion at &-4.84^, for an inviscid Euler 
calculation. This Euler solution uses a 121x35 C-grid for the computation. 
Similar results were produced for ensuing angles of attack. Figures 35 through 
39 show pressure distributions at a=4.84°, a@=6.72°, a@=10.80°, @=12.83°, and 
0=15.54°. Surprisingly, the inviscid code predicted the suction peak more 
accurately than the N-S code for all angles of attack. Several explanations are 
offered to account for this unexpected result. First the 161x65 C-grid may not 
be fine enough to predict the magnitude of the suction peak and the N-S code 
may have too much numerical viscosity, thus in effect performing calculation at 
a lower Reynolds number. The flow at the leading edge may in fact be laminar, 
while the computed solution assumed turbulent flow everywhere. Figure 40 


shows the Mach contour and Figure 41 shows the density contour at a=10.8". 
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Figure 34. Pressure Distribution (Inviscid), Ramp Motion, M=0.3, 
(24.87, kz0.01272, Ке-2.7х106 
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A detail of the flow field generated by Ns2.f at &-12.83' 1s shown 
in Figure 42. The flow appears to be smooth and attached (the dark line shows 
the boundary layer edge). However, magnification, as shown in Figure 43, 
reveals the start of a slight reverse flow region at approximately 10% chord. At 
a=15.54°, Figure 44, the boundary layer is much thicker with reverse flow 
profiles from approximately 10% chord to the trailing edge. 

Figures 45 and 46 present the computed C] vs « and Cg vs & 
curves, respectively. The Navier-Stokes code underpredicts lift and pitching 
moment coefficients whereas the inviscid panel method gives significantly 
better lift predictions. However, it fails to predict the pitching moment 
coefficient for angles greater than 10°. A more comprehensive investigation of 


the dynamics of a rapidly pitching airfoil can be found in Grohsmeyer [Ref. 19]. 
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3. Harmonic Airfoil Oscillation Using the Baldwin-Lomax 
Turbulence Model 
Next, the code was applied to sinusoidal airfoil oscillations about the 
quarter chord point. These time periodic solutions were obtained for the second 
cycle because the results obtained for the third cycle were indistinguishable 


from ones of the second cycle. Figure 47 shows the computed angle of attack 


history. 


Angle of Attack 


Steady State First Cycle Second Cycle 





0 10000 20000 


Figure 47. Harmonically Oscillating Airfoil 


Here, AQ is the mean angle of attack, and Aj is the amplitude. For 
each oscillatory case a steady-state solution was obtained for the minimum 
angle of attack during the cycle. The iteration counter was set to zero as the 
sinusoidal motion was time shifted half a cycle so that the start of the motion 
begins at the minimum angle of attack. This ensures a smooth transition for the 


converged steady state solution to the start of the oscillatory motion. Ten 
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thousand (10,000) iterations were computed per cycle, with a timestep equal to 
approximately 0.05 to 0.01. 
a. Case 6. Q(t)=Ag°tA]°sin(@t) 

All the computations presented in this section are computed with 
Ns2.f using the Baldwin-Lomax turbulence model. Figures 48 and 49 show the 
effect of Mach number on the lift and moment hysteresis loops for small 
amplitude oscillations about a mean angle of 12 degrees with an amplitude of 
3.0 degrees, at a Reynolds number of four million and at a reduced frequency 
of 0.1. Note that the hysteresis loops for Moo=0.4 oscillate wildly and third 
cycle computations do not coincide with results from previous cycles. For the 
M.o=0.3 computation the moment hysteresis loop indicates a stable condition. 
The effect of amplitude is shown in Figures 50 and 51. Figures 52 and 53 also 
show the effect of amplitude for a larger mean angle of attack, Ag=14.0", at the 
same Mach number, reduced frequency, and Reynolds number, as before. Both 
solutions oscillate in a manner that is not repeatable if the motion is allowed to 
continue for additional cycles. For this comparison the maximum angle of 
attack attained during a cycle exceeded 16.5". A comparison of the first and 
second cycle lift and moment coefficient loops is shown in Figures 54 and 55. 
Note that the second cycle lift curve starts at a slightly greater lift coefficient 
than in the first cycle. Aside from this anomaly, the loops are in good 
agreement but exhibit no instability. 

Figure 56 shows the computed lift and moment loops for small 
amplitude oscillations about a mean angle of 13 degrees with an amplitude of 


2.5 degrees, at a Reynolds number of two million, for three values of reduced 


9] 


frequency. It can be seen that all the moment loops computed in this section 


using the Baldwin-Lomax turbulence model indicate torsional stability. 
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Figure 48. Effect of Mach Number on Unsteady Motion, 
a(t)-12^x3'sin(OGt), k-0.10, Ке-4х106 
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Figure 49. Effect of Mach Number on Unsteady Motion, 
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Figure 50. Effect of Amplitude on Unsteady Motion, M=0.3, 
.(1)-12:3А1 Ssin(@t), k=0.10, Re-4x106 
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Figure 55. Effect of 1S and 2nd Cycle on Unsteady Motion, M-0.3, 
Q(t)=13+2°sin(wt), К-0.10, Ке-4х106 
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b. Case 7. a(t)z13^42.5^sin( ою), К-0.2, Ке-4х106 

In an effort to determine if the flow field was in fact predicting a 
trailing edge vortex that would shed during the oscillation cycle as seen in the 
experimental results, instantaneous particle traces and the corresponding 
velocity fields were plotted for the downstroke at a Mach number of 0.3, a 
mean angle of attack of 13 degrees, an amplitude of 2.5 degrees, a reduced 
frequency of 0.1 and a Reynolds number of four million. Figures 57 through 61 
show the downstroke portion of the cycle beginning at @=15.3° and ending at 
a=14.5°. A recirculatory region is seen in Figure 57 which grows slightly, as 


seen in Figure 58. This region then appears to diminish in size and intensity. 
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c. Case 8. a(t)213^42.5^"sin( Qt), k-0.1, Re-2x106 

The plots generated for Case 7 did show a recirculatory region, 
but a full understanding of its development and structure could not be gained 
from only a portion of the oscillatory cycle. Therefore in Figures 62 through 80 
a complete cycle is shown for the same flow condition as in Case 7 with only 
the reduced frequency being increased from 0.1 to 0.2 and Reynolds Number 
decreased from four million to two million. A recirculatory region is seen to 
form at a=14.9° on the upstroke, which continues to grow and intensify. This 
recirculatory region reaches its largest value at @=15.3° on the downstroke and 
then diminishes and finally the flow becomes completely reattached. It is 
important to note that even for this relatively fast oscillation the recirculatory 


region did not shed into the wake. 
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Figure 72. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.3 IL 
a(t)=13°42.5°sin(@t), k=0.10, Re=2x10 


119 





NORMAL | ZER iria T 108. PRESSURE -e 


— — 
— — = — — — = . 8213.0, ñz2.5. кі0.1, B-L м 410оип51 гоке ) 
a RET "MEC P n e 
= Өг сс чие s — => 2 По == 
CONTOUR LEVELS ES les S 0.300 nat 





Figure 73. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.1° (Dovnstroke 
«(t)=13°+2.5°sin(®t), К-0.10, Ке-2х10 


120 






=== =РВЕ ТЈ - 0 = 
Е 


i ms 


NORNAL I ZER ÜÓEWIAGT IOS PRESSURE m 


_ 9213,0, 822.5, к=0.1, B-L tet -100unst roke ) 7 


w 
Å- тыь 
“ж. 


~ 
в“ 
E Л -— 
— = < - س‎ Te -— 
И 


“CONTOUR LEVELS в ИН ~ 
(1. ang ——— _ 


~ 
m € 0.300 HRCN-— 
ےر‎ 97 — =  — == 


xe ~ !4.90EG ALPHA 
л. k “эрх 104686 
т" и rcl 8 жы ы — 


ИИ = 1. 69x tas? TINE 
E mr RT: И DE a 


| " e 1644 т 
E cl ха 





L, 


Figure 74. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.9 етос. 
a(t)213 22.5 sin(Ot), k20.10, Rez2x10 


124 





NORMAL | Zea GEWWART ION PRESSURE - 


~. E 
ше e æ e ص‎ 9=13.0, 822.5, кгр. l, B-L t7 -—Qounst roke ) 
acp LL а — -— з. E 
a ee “ы т. MELLE 

-CÜNTOUR LEVELS SAOR 2277 0.300 RC — 

pp АИ NS =. 14. 7066 ALPHA 

"СТО pree -æ "22x 10**6Re 

= LM e. ~ -.. m 16464 CR 
S Беча ар c = = 





Figure 75. Instantaneous Particle Trace and Velocity Field, 
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Figure 77. Instantaneous Particle Trace and Velocity Field, 
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4. Effect of Turbulence Modeling 
a. Case 9. a(t)=9°+5°sin(@t), k=0.2, Re=4x106 

Figures 81 and 82 show a comparison of the lift and moment 
hysteresis loops using the Baldwin-Lomax turbulence model for an oscillation of 
5 degrees amplitude about a mean angle of attack of 9.0 degrees with 
McCroskey's experimental data for a Mach number of 0.3, a Reynolds number 
of four million, and a reduced frequency К=0.2. It is seen that there is a 
substantial difference between the computed and the experimental lift 
hysteresis loops. This difference becomes even more pronounced when the 
computed and experimental moment hysteresis loops are compared with each 
other. In an effort to understand the failure of Ns2.f to predict the experimental 
pitching moment hysteresis loop, the pressure distributions for several upstroke 
and downstroke angles of attack were compared to measured data from 
Reference 1. Figures 83 through 87 show pressure distributions for a@=11.9" оп 
the upstroke portion of the cycle; and @=13.7°, о-13.07, a=10.5°, and a=8.9° on 
the down stroke cycle. U2diif.f overpredicts the suction peak at a=11.9° on the 
upstroke and Ns2.f underpredicts the peak. Even though the values are off, the 
general shape of the distribution is well predicted. However, on the downstroke 
at @=13.7", a significant difference between the computed and the measured 
distributions is observed. The measured pressure distribution shows a lower 
suction near the leading edge but higher suction near the trailing edge. Neither 
of these changes are predicted by Ns2.f. As the oscillation continues on the 
downstroke the measured distribution tends to flatten out creating a "plateau". 
Just after &—8.9' on the downstroke, the flow begins to reattach and the positive 


pressure difference at the trailing edge becomes smaller. The failure of Ns2.f to 
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predict the experimental loop results from the inability of Ns2.f to accurately 
compute the pressure distributions. Since the pitching moment is an integrated 
quantity, the size of the area under the pressure distribution aft of the quarter 
chord must be greater than the area forward of the quarter chord to compute a 
negative or downward pitching moment. Ns2.f clearly does not predict the 


experimental pressure distributions during the downstroke. 
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b. Case 10. a(t)=10°+5.5°sin(@t), k=0.1, Re=4x106 

The measured loop consists of two subloops, a clockwise 
“unstable” subloop and a counterclockwise “stable” subloop. It is seen that the 
Navier-Stokes calculations in combination with the Baldwin-Lomax turbulence 
model fails to predict the destabilizing clockwise pitching moment loop. For this 
reason the sensitivity of the computed loops to different turbulence models was 
studied next. Figures 88 and 89 show the computed lift and pitching moment 
loops for oscillation about a slightly higher mean angle of 10 degrees at an 
amplitude of 5.5 degrees. The Mach and Reynolds numbers are again 0.3 and 
four million, and k was decreased to 0.1 to minimize the effects of reduced 
frequency. It can be seen that the Baldwin-Lomax, Johnson-King, and the 
RNG turbulence models produce significantly different hysteresis loops. The 
RNG model produces relatively good agreement with the measured lift 


hysteresis loop but fails again to predict the destabilizing moment subloop. 
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c. Case 11. a(t)211^X5'sin(Qt), k=0.1, Re=4x106 

A solution was also computed for oscillations about a mean angle 
of 11 degrees at an amplitude of 5 degrees, the other parameters left 
unchanged. Figures 90 and 91 show the computed lift and moment hysteresis 
loops. As before, the calculations fail to predict the expected destabilizing 
moment subloop. The reason for this failure can be better appreciated by 
visualizing the separated flow structures computed with these turbulence 
models. Figures 92 and 93 show the computed flowfield with particle traces 
and the velocity field. Figure 92 corresponds to the instant when the airfoil 
oscillates through an incidence of 15 degrees during the downstroke. Different 
turbulence models produce substantially different recirculatory flow patterns on 
the upper surface near the trailing edge. The Baldwin-Lomax model produced 
the smallest recirculatory region while the RNG model produced the largest. In 
addition, Figure 93 (az15^) shows the relative magnitudes of the velocity 
vectors. Note, near the surface at the trailing edge, the magnitude of the 
velocity vectors in the reverse flow region calculated by the Baldwin-Lomax 
model are much smaller than the region calculated by the Johnson-King model 
or RNG model. As the cycle continues, Figures 94 and 95 (a=14") show the 
recirculation region to shrink slightly and loose intensity for the Baldwin-Lomax 
model, while the regions as calculated by the Johnson-King and RNG continue 
to expand and intensify. As noted in Case 8 the recirculatory flow region 
grows and then vanishes again during the course of the oscillation. Even 
though the structure at the trailing edge resembles a vortex it does not shed 
from the trailing edge. The occurrence of the destabilizing moment loop 


appears to be closely associated with the vortex shedding from the trailing 
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edge, an event which occurs soon after the static angle is substantially 
exceeded during part of the cycle. Although the static stall angle is significantly 
exceeded in the present calculations the numerical procedure along with the 


turbulence modeling used in the calculations appears to be unable to produce 


the anticipated shedding process. 
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EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
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Figure 92. Instantaneous Particle Trace, Oscillatory Motion, M=0.3, 
a=15° (Downstroke), a(t)=11°+5’sin(wt), k=0.10, Re=4x106 
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EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
TRAILING EDGE FLOWFIELD 
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Figure 94. Instantaneous Particle Trace, Oscillatory Motion, Mz0.3, 
a=14° (Downstroke), a(t)=11°+5°sin(@t), k=0.10, Re=4x106 
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V. CONCLUSIONS AND RECOMMENDATIONS 


In the preceding investigation a computationally efficient, fully factorized, 
two-dimensional Navier-Stokes flow solver was utilized to predict steady and 
unsteady flow solutions about a NACA 0012 airfoil. Comparisons between a 
steady viscous/inviscid interaction method code using the Cebeci-Smith 
turbulence model and the Navier-Stokes code using the Baldwin-Lomax 
turbulence model show close agreement up to just prior to the static stall angle 
of attack. After this point, the viscous/inviscid interaction method code fails to 
accurately model the flow separation. 

When applied to the unsteady problem of airfoil stall flutter in compressible 
flow, the Navier-Stokes code shows that the modeling of the recirculatory flow 
region with current turbulence models fails to capture the essential physics 
which governs the onset of stall flutter. Comparisons of the computed results 
with available experimental data indicates that even though the lift response is 
fairly well predicted, the computation of the pitching moment hysteresis loops is 
very sensitive to turbulence modeling. Results computed with several current 
models are in good agreement whenever the steady stall angle is exceeded only 
slightly. However, they fail to capture a vortex shedding process that may 
contribute to the onset of stall flutter. 

Therefore, further detailed studies of improved numerical schemes and 
turbulence models as well as viscous/inviscid interaction approaches are 
required to improve the prediction of the unsteady flow separation and vortex 
shedding phenomenon. Further computational work with the full Navier-Stokes 


equations instead of applying the thin-layer approximation is recommended. In 
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addition, local grid refinement studies that focus on the the leading and trailing 
edge upper surface may lead to the accurate prediction of the extremely 


complex vortex development and shedding process. 
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APPENDIX A - FLUX JACOBIAN MATRIX 


The Jacobian matrices A = ЈЕ / од апа В-Ф/ O are given by Aor В- 
КОГ + К1А + k2B, where A = ðE/ðq апа В = OF/dq are the usual Jacobian 


matrices of the Cartesian flux vectors. The Aor B matrix is 


gu | 
б -(у-1) ғ UE | 
| tk заг Киа: Еш : | x 

- OF 

ы SVK U +k,v) Кум -(y - 2)k,v * | ТШ | 
| tke AY- Dk; k*kuckv — 7 77 | 
| | 
| (ku +k,v) * ly га. Iy(e/ p) - e? ]k, - y (ku + kv) | 
[-76/2)+29] (у-1\(ки+к.ўн (у-1кш+кд)у +, | 


where 62 = 0.5(y - 1)(u2 + v2), kg=Et, k1-5Ex, k2-£z for A, алко С kK =x. 
k2=űz for В. 
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APPENDIX B - PITCHING MOMENT DERIVATION 


Figure 96 shows the nomenclature for the integration of the pressure and 
shear stress distributions over a two dimensional body surface as generated by 
the Navier-Stokes code Ns2.f. For unsteady motion the airfoil is rotated with 
respect to a fixed inertial frame of reference (x,z: Laboratory frame of 
Reference). Note, for unsteady motion, as the airfoil rotates through an angle of 
attack, all grid G=1...Imax, k=1...Kmax) are rotated about the desired pivot point 
with respect to the fixed x,z coordinate system. As the solution is advanced in 


time the flow quantities for this new grid location are updated for each point. 


Е upper 1131 


A 
. . | 
Az-z(i*1) - 20) ~ T.E. lower, i=31 
k • • 2 


Ax=x(i+1) - x(i) 


Z - coordinate 


x - coordinate 


Figure 96. Pitching Moment Derivation 
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PITCHING MOMENT DERIVATION GENERATED BY MATHCAD 


1. Read tabulated data generated from Ns2.f. Data computed at a=0° and a Re=4x10**6. 


Data := READPRN (up) Data : = READPRN (low) i-um ll 5 
u 1 
<0> <0> 
х : = Data х : = Data 
u u 1 1 
<> «1» 
2 := Data 2 >= Data 
u u 1 1 
<2> <2> 
ср :ж раба ЕР := Data 
u u 1 1 
<3> <3> 
T := Data Л := Data 
u u 1 


2. Compute Ax and Az for upper and lower surface. 


Ах бах - x Ax. “= ин - x.l 
u u u 
i 1 с: i i (1-1) 


Ат =? - 2 Ат.1 = Z2] = z.1 


1 1. 1-1 1 1 (2-1) 
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ı Computed the x and z moment arms for the upper and lower surface. 


хх „= 22 > = 


| Computepitching moment coefficient about quarter chord by applying equation 1.11 from [Ref. 
lp. 17]. 


зер °cos |8 кт "sin|8 "XX 

u u u u u 
қ т 1 1 I Qc cuis 
u 


+ т *"sin|8 “хх 
L 1 1 
1 1. li QE LENS 


. The computed pitching moment coefficient from Ns2d.f (Computed on a Cray Super-Computer) 
'а5 Ст=0.0006, as compared to the calculations made on a Macintosh SE/30 where Cm=0.005. 
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Upper Surface Data: 


X 


m. 
=2. 
2r 
=. 
r: 


202 
= 
=2 
=2 
E 
21 
= 


-4 


~ ~] ДА ОУ бу O O (л S +> > а (ЈОЈ ОЈ ММ ММ нь -1. tS) b 


5е-1 

497e-1 
479е-1 
447е-1 
404е-1 


.346e-1 
.273е-1 
.187е-1 
.088е-1 
.975е-1 
ВЕ 
.713е-1 
-1. 
i 
-1, 
-1. 
— 
-6. 
. 38e-2 
-2. 
236-1 
.5е-2 

. 95е-2 
.45е-2 
20227651 
.264e-1 
e 
.8е-1 
.073е-1 
.348е-1 
.626е-1 
.904е-1 
.182е-1 
.461е-1 
Te 
.014е-1 
.287е-1 
.557е-1 
.823е-1 
.085e-1 
.342e-1 
.593е-1 
.838е-1 
.076е-1 
. 306e-1 
526651 
.742e-1 
.946e-1 
.141е-1 
.325е-1 
.499е-1 


564e-1 
403e-1 
23e-1 
047e-1 
54e-2 
yc 2 


17e-2 


7. 
zu 
3 


ON -J F2 P2 F2 F2 о ОЈ CO CO uS d» uS oS uS OI UI O1 Qi Cn Q Q Q OT tr л Q Q Q Ui On Cn uus uS Q) CO CO RO МЮ NO ES F9 «J 


-C Tau 
. бе- 3 -1.0319е+0 Зе-3 
е-3 -9.841е-1 -7.4е-3 
.бе-3 -6./98e-1 JI EL 
.2е-2 -4.091е-1 -1,64ё-2 
.65е-2 -1.766e-1 -17 54e-2 
.09е-2 9.93е-2 -1.19е-2 
.49е-2 2.286e-1 -8.4e-3 
.88е-2 2.769e-1 -6.3e-3 
.26e-2 3.316e-1 -4.9e-3 
.61е-2 3.715е-1 -3.8e-3 
. 94e-2 3.951e-1 -3e-3 
.26e-2 4.105e-1 -2.4е-3 
.55e-2 4.199e-1 -1.9e-3 
:81е- 7 4.245e-1 -1.5e-3 
.05е-2 4.248e-1 -1.2e-3 
.26e-2 4,221е-1 -1е-3 
.44e-2 4.168e-1 -7е-4 
. бе-2 4.094е-1 -бе-4 
.73е-2 3.996e-1 -4e-4 
. 82e-2 3.904e-1 -Зе-4 
. 89e-2 3.778e-1 -2е-4 
. 93e-2 3.628e-1 -1e-4 
. 94e-2 3.494e-1 безо 
. 92e-2 3.343e-1 1е-4 
Ec о 3.184e-1 1е-4 
. 81е-2 3.031е-1 2е-4 
.72e-2 2.875e-1 2e-4 
. бе-2 2.71е-1 2е-4 
.46e-2 2.544e-1 3e-4 
. 36-2 2.378е-1 3e-4 
13672 22116-1 3e-4 
.93е-2 2.045е-1 Зе-4 
.72е-2 1.884е-1 Зе-4 
. 5e-2 1.716e-1 4e-4 
.27e-2 1.559e-1 4e-4 
20276—7 1.421е-1 4е-4 
.77e-2 1.252e-1 4e-4 
.5е-2 1.067е-1 4e-4 
.23е-2 9.05е-2 4e-4 
. 96e-2 7.46e-2 4e-4 
.68е-2 5.68е-2 4e-4 
.4e-2 3.86e-2 4e-4 
12652 2.06е-2 4e-4 
.84е-2 1.3е-3 Зе-4 
.56е-2 mU S 2 3e-4 
. 29e-2 -4.42e-2 3e-4 
.01е-2 -7.02е-2 Зе-4 
.4e-3 -9.2e-2 2e-4 
· је=3 -1.441е-1 1е-4 
.2е-3 -1.088е-1 0е-0 
е+0 -2.34е-1 Ое+0 
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Lower Surface Data: 


X 


-2. 
-2. 
.458е-1 
.419e-1 
.366e-1 
2 26-1 
.215е-1 
.12е-1 
.012е-1 
.891е-1 
не! 
.613e-1 
.456e-1 
.288e-1 
.109e-1 
.19е-2 
.2е-2 
.11е-2 
. 94е-2 
.8е-3 
,66e-2 
.06е-2 
3озе-о 
.06е-2 
.164е-1 
.427e-1 
.694е-1 
.963е-1 
.23бе-1 
яоТте- 1 
. 785е- 1 
.061е-1 
.336e-1 
.61e-1 
.883e-1 
ul54de-l 
.421e-1 
.685e-1 
.945е-1 
.199е-1 
.448е-1 
.69е-1 
.926е-1 
«І Зще- 1. 
4 /5е=1 
noo 7 e— 1. 
.789e-1 
. 982e-1 
.165е-1 
.338е-1 
.499е-1 


J لہ‎ -J Gy OY. O03 0) 60Y UU! л сл > из ы ыы, 69 CO CO CO MON NF Pe Po oy > i 


5e-1 
486e-1 


Z 


-1l. 
-6. 
EU. 
-1. 
-1. 
-2. 
-2. 
- 3. 
-3. 
-3. 
.16е-2 
-4, 
-4, 
-4, 
-5. 
-5. 
-5. 
E 
-5. 
— 
-5, 
-5. 
-5; 
-5. 
-5. 
-5. 
zi 
-5. 
-5. 
-5. 
=e 
-4, 
-4, 
-4, 
-4, 
-3. 
-3. 
нэ 
-3. 
-2. 
-2. 
-2. 
-2. 
EL. 


-4 


T 
zi 


6e-3 

2е-3 

06e-2 
51e-2 
96e-2 
37e-2 
76e-2 
14e-2 
5e-2 

84e-2 


45е-2 
136-2 
97e-2 
19e-2 
39e-2 
55e-2 
69e-2 
79e-2 
87e-2 
92e-2 
94е-2 
93е-2 
9е-2 

84e-2 
76e-2 
65e-2 
52e-2 
37e-2 
2е-2 

02е-2 
82е-2 
бе-2 

зоба 
14е-2 
89е-2 
64е-2 
37е-2 
11е=2 
84e-2 
57e-2 
29e-2 
02e-2 
75e-2 


.48e-2 
.21e-2 
-9. 
-6. 
— 
-2. 

Qe+0 


5e-3 
9e-3 
4e-3 
le-3 


| О OY CO XO НЗ HP P) P) 2 i Му Му М М NO Ñ) Q) O) O) Q) CO CO CO ць из из PB PW WWD NN 


3e-3 
69e-2 
06e-2 
47e-2 
64e-2 
491e-1 
101e-1 


.0319е+0 
.832е-1 
.815е-1 
.554е-1 
. 4е-2 

.063e-1 
.624е-1 
.142е-1 
.619е-1 
.892e-1 
.067е-1 
.178е-1 
.239е-1 
1256621 
.237е-1 
.191e-1 
.126e-1 
. 035е-1 
.932е-1 
.84е-1 

.681е-1 
.527е-1 
.419e-1 
.249e-1 
.084е-1 
.941е-1 
.777е-1 
.612е-1 
.448е-1 
.283е-1 
„өш il 
Tee l 
.794е-1 
.631е-1 
. 49е-1 

.342е-1 
.16е-1 

.88е-2 
.37е-2 
.736-2 
.95е-2 
.18е-2 
.43е-2 
E 
-2. 
шог 
-7, 
-9, 
—1. 
-1. 
-27146-1 


Ое+0 
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Таџ 


3e-3 


4. 
2: 


9e-3 
3e-3 


-2е-4 


-2. 
-2. 
— 
=2, 
22. 
—1. 
-1. 
zu. 
du 


2e-3 
8e-3 
7e-3 
5e-3 
2e-3 
9e-3 
7e-3 
4e-3 
2e-3 


-1e-3 
-8e-4 
-7е-4 
-5е-4 
—4е-4 
– Зе-4 
-2e-4 
-16-4 
0е+0 
1е-4 
1е-4 
2e-4 
2e-4 
Зе-4 
Зе-4 
3e-4 
4e-4 
4e-4 
4e-4 
4e-4 
4e-4 
4e-4 
4e-4 
4е-4 
5e-4 
5e-4 
5e-4 
5e-4 
4e-4 
4e-4 
4e-4 
4e-4 
4е-4 
4e-4 
3e-4 
2е-4 
1е-4 


APPENDIX C - USING THE PROGRAMS 


A. DIRECTIONS FOR THE EXECUTION OF NS2.F 
1. Steady State Case 


In order to generate a Cp, C], or Cm plot from Ns2.f there are many steps, using 
six different programs that must be taken: (1.) grape (grid generation program), (2.) rotgr 
(grid rotation program), (3.) Ns2.f (Navier-Stokes code), (4.)vi (editor), (5.) plot3d 
(flow visualization program), and (6.)xyplot (a simple x-y plot program). 

A.) Ns2.f is a fortran code and will take hours to run even on the Stardent mini- 
super computer. To run this code, an executable file (ns2), an input file ns2.in (actual 
name of input file is users choice) and a grid file FORO11.DAT (Stardent) or fort.11 (on 
other computers) must be in the directory. 

B.) To obtain a grid, a grid generation program such as grape must first be 
used. Grape output places the leading edge of the grid at the reference axis at 0.0° AOA. 
Now, the grid must be rotated to the desired AOA by using the program rotgr. The input 
file (output from grape) for rotgr must be FOROO1.DAT and the output is FORO02.DAT. 
The output from rotgr (FOROO2.DAT) must be renamed to FOROI 1.DAT for use in ns2. 

C.) Next edit the input file to reflect the desired flow conditions. Any editor will 
do, however vi is perhaps the most common editor and it is found on most unix based 
machines. Enter Mach No. and AOA, check smoothing factors, set time step, Courant No. 
and initial no. of iterations (small value- for first run). Check restart false for first run and 
make sure No. point of your grid match input values. 

D.) The code is executed as follows: ns2 <ns2.in> ns2.out 

The first ten iterations should take 1-10 minutes, depending on the type of 
computer used. Several output files will be generated: 


TABLE 3. NS2.F OUPUT 


COMMENTS 
convergence data 
Ci da 

residuals 

Ci, Ca, Cn data 
Cp data (io be used with xyplon 
FORO31.DAT: Flow data (to be used with plot3d) 






















E.) Re-run the code for a greater number of iterations to achieve a 
converged solution. When the change in the max residual (FOROO7.DAT) has dropped 


1072, the solution has converged. Ensure restart is set to true. The CFL can also be 
adjusted to obtain a converged solution more rapidly. 


158 


МАСН, 
02:00, 
EDAX, 
nono, 
DT, 
0-005020. 
КОЈЕГ. 
те, 
TIMEAC, 
1; 
Vasc, 
true, 
DFEL, 
al 
ІА1, 


ШӘШШӘЕШ ЕП; 1550,157071590, 1610,1630,1630, 


false 


UNSTST (If true Time - 
steady state, 


a. Sample Input to Ns2.f 
ALFAO, АТЕА1, ALFARE, REDFRE, REYNOLDS 
0:00: 20-10) (200. 0700; D 
Брате EDAX, ED4Y, ED 
GOO OOO, 020, 200 
COURT NITER FE NEWTIT 
21009 307 2 
OS CIE; RAMP, NPER, ШЭНТЕТ 
false, false, 10000; _ = 
IMPLBC, EXPLBC, СВЕ 
false, true, false 
BLTM, JKIM, RNGTM 
true, false, false 
ТТЕП) 
ШЕЛІ 
БАЛА ТАЯ TES ЖАС, ТА), ТАЗ, ІАЗ, ТА1О 


for proper recording of unsteady motion) 


Mach 
AlfaO 
Alfal 
Alfare 
Redfre 


Reynolds 


ED2X 
ED 
ЕР4х 
ED4z 


ED 
DSPEC 


Dt 
Cour 
Niter 
Newtit 


ВЕЕТ 
OS EL 
RAMP 
NPER 
TONPET 


ТЇМЕАС 
IMPDBC 
EXPLBC 


Free stream Mach number 


1805 


Angle of attack, also mean angle of attack for 


Amplitude of Oscillatory motion 
Reference angle 

Reduced frequency k = omega * c 
Reynolds number Re = cU/n 
explicit 


X-direction 2nd order 


Z-direction 2nd order explicit 


се п лла order explicit 


Z-direction 4nd order explicit 
Scaling of Implicit smoothing 


Spectral radious parameter 


стер 

tonrcantenumber Cu — dt * L max 
Number of Iteration in this run 
Newton subiteration within each 


Restart 


20 


smoothing 
0525 
smoothing 
827110, 
smoothing 


smoothing 


timestep 


( 
< 
( 
< 
( 


( 


е2х 
е2х 
е22 
e2z 
е4х 
e4x 
e4z 


0.,Set TRUE only when unsteady motion starts from steady 
TRUE for steady-state ok, but for unsteady restarts must be FALSE 


unsteady 


[ I та | 


| 
сае ос 


ейт = 


ЕЕ, ПОТОП А (Е) = AO + Al * sin ( k * M * t ) 


Ramp motion 


Number of time steps in one period of oscillation, 


Time accurate Тасс=1 апа for Jacobian Scaled Dt, 


Implicit wall BC Treatment 
Explicit wall BC Treatment 
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subsonic, 
transonic) 
subsonic, 
transonic) 
subsonic, 
transonic) 
subsonic, 
transonic) 


dt-T/Nper 
Time shift in radians to start oscilation for any a(t) 


Тасс-0 


УТС : Only laminar viscosity 


TURBL : Baldwin-Lomax eddy viscosity 

JKTM : Johnson-King eddy viscosity 

RNGTM : RNG eddy viscosity 

Т ГЕТ ITEU : Lower and Upper trailing edge I locations 

IAl etc. : Write out (IAn/100) angles of attack in units 
З 7 


Read grid from unit ғогЕ.11 апа the flow шт 1: 0201 
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b. Sample Output from Ns2.f 
Grid dimensions : 161x64 


Mach = 0. 30000 
AO = 0.00000 Al = 0.00000 k = 0452722 
Ке zl 220800ЕТ-097 Dt - 0.00200 mm —7009)0900000 
Dimpl = 2.00000 
р2х = 0.00000 D2z = 0.00000 
bax = (204:000” 2727 = 0.03000 
Restart = ЈЕ 
Oscil = F 
Ramp = В 
Itel - 31 Iteu = 131 Ile = 81 
Timeac = T 
BC impl = Е 
BC expli - T 
Cour NONO 
L max = 232096.4825513 
Dt = 0.00904796 
ТЕ drmax dumax dvmax demax qe kr 
Богт ш 71716Е-О?” (1: :563395Е-07 0.129865Е-06 ДЕДІ? 63 
Ni— J А 05%0:1711341E-08 02803250Е”08 0.219128Е-07 
Alfa(t) = 0.000 C1 = 0 007 У ве са = 0.003260 Cm » -0.000692 
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2. Unsteady Case 

The execution of an unsteady case is essentially the same as for the 
steady-state case. 

A. First, a well converged steady-state solution must be obtained. 
This solution then becomes the starting point for the unsteady solution. 

B. The restart is left set to true, however, to reset the time counter 
to time=0, the UNSTST variable must be set to true for several iterations and 
then back to false for normal time counting. 

C. As the unsteady motion begins to march in time Flow data and 
Grid position data may be stored at user determined angles of attack. To store 
flow data at @=12.5°, a@=13.5°, a=14.5°, etc... enter the data shown in Table 4. 
into the input file. 





D. The Flow data and the Grid data at each œ will be stored in data 
files FORO21.DAT through FORO30.DAT. The Grid and Flow data can be 
separated using a simple fortran code that reads the data and then writes the 
data to two separate files. Splgf.f inputs the desired file ID and outputs a Grid 
file and Flow file. 
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МАСН, 
02900, 
ЕР2Х, 
O00; 

ID 
020220 
RSTRI, 
true, 
TIMEAC, 
1, 
MISC. 
true, 
ITEL, 
E 
ТАТ, 


ШЕГІ ФИО ood, 570, 1590, 1610,1630,1630, 


false 


a. Sample Input to Ns2.f 


ALFAO,  ALFAl,  ALFARE, КЕРЕКЕ, REYNOLDS 
0. 2 15254. 05000701272, 2. О 
RIDE EDAX, ED4Y, ED 
(ОД и а) (27205508 210 
COURS МИРЕ А НИ ЕЙИГТТ 
2100, SÛ, al 

CSCI, RAMP, NPER, TSHIFT 

false, true, 30007 -0.5 

IMPLBC, EXPLBC, CIRCOR 

false, true: false 
BLTM, JKTM, RNGTM 
true, false, false 
ТТЕП 

151 
ОА, Па TE TAG TA? ТАЄ, ТАЗ, ІА10 


ew 


UNSTST (If true Time = 0.,Set TRUE only when unsteady motion starts from steady 
steady state, TRUE for steady-state ok, but for unsteady restarts must be FALSE 
for proper recording of unsteady motion) 


Mach 
А1Ға0 
Alfal 
Alfare 
Redfre 


Reynolds 


ED2x 


ED2Z 


EDAX 


17 
ЕР 
ISPEC 


Du 
Cour 
Niter 
Newtit 


вт 
OSCTL 
RAMP 
NPER 
TSHTET 


TIMEAC 
IMPLBC 
EXPLBC 


Free stream Mach number 
Angle of attack, 
Amplitude of Oscillatory motion 

KETE NCE angie 

Reduced frequency k - omege * c / 2U 
Reynolds number Re = cU/n 


smoothing 
25 
smoothing 
O10 
smoothing 


X-direction 2nd order explicit 


Z-direction 2nd order explicit 


X-direction 4nd order explicit 
Z-direction 4nd order explicit 
Seo OE IMOIICIL smooctbasng 


Spectral radious parameter 


smoothing 


Time step 

Courant Бег Си = ЧЕ ~ пах 

Number of Iteration in this run 

Newton subiteration within each timestep 


Restart 
Oscillatory motion A(t) - 
Ramp motion 


also mean angle of attack for 


e2x 
e2x 
e2z 
e2z 
е4х 
е4х 


unsteady 


ЕЛ ШЕЛ II 


Qa C ea] 


е42 = 


Number of time steps in one period of oscillation, 


шиг КМ * Lt ) 


subsonic, 
transonic) 
subsonic, 
transonic) 
subsonic, 
transonzc) 
subsonic, 


dt=T/Nper 
Time shift in radians to start oscilation for any a(t) 


Time accurate Tacc=1 and for Jacobian Scaled Dt, Tacc=0 


Implicit wall BC Treatment 
Explicit wall BC Treatment 
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VES : Only laminar viscosity 


TURBL : Baldwin-Lomax eddy viscosity 

JKTM : Johnson-King eddy viscosity 

RNGTM : RNG eddy viscosity 

MTEL, ITEU : Lower and Upper trailing edge I locations 

IAl etc. : Write out (IAn/100) angles of attack in units 
P t 


Read grid from unit fortil and he fle. from опа 00-70 
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b. Sample Output from Ns2.f 


Grid dimensions : 161x64 


Mach -з 0.30000 
Ao - 0.00000 Al = 16.00000 k = 021212 
Ке = 0.2700Е+07 Dt = 0.01829 Cu= 15.00000 
Dimpl = 2.00000 
Dyk = 0.00000 D2z = 0. 00000 
р4х = 0.03000 D4z = 0.03000 
Restart = Tr 
OSE = Е 
Катр = T 
Itel = 31 Iteu = 131 Ile = 81 
Timeac = 1 
Вор = Е 
BC expl - T 
Cour = А209 7352745382 
E MAS Е = РАО озо. 9275225 
Dt = 0.01829000 
It drmax dumax dvmax demax ir 
шэг 02249290Е“06 0:1021934Е-05 0.190543Е-05 0.174200Е-05 1219 
Ni= 1 0.113997E-06 0.302620Е-06 0.448250Е-06 0.196797Е-06 


ЕО АЕ О omega = 0.007632 phase = 0.043322 
1.12069 Са - 0 157214 Cm = 0.003384 


ог: 202107038Е4-05220-191827Е-05 0.158330Е-05 TTS 
Ni= TID TEE OG 0-297453E-06.320.450476E-06 0.197668Е-06 


it = 2000 Time = 36.580000 omega = 0.007632 phase = 0.044433 
Alfa(t) ج‎ ш  32217- 1.455050 Cada = 0.146989 Cm - 00015295 
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КЕ 
63 


63 


B. DIRECTIONS FOR THE EXECUTION OF U2DIIF.F 


]. Unsteady Case 


A. A complete description of the input and output parameters can be 


found in Teng. 


TABLE 5. U2DIIF.F INPUTS/OUTPUTS 


FILE COMMENTS 


FOROOLDAT 
Оза ош 


B.) For a unix based machine the code is executed as follows: 





U2diif 2». U2diif.out 


Note: The original Fortran code can be modified to store desired output in 
separate data files by inserting a WRITE (9,*) statement inside the required do- 
loop. (Example: Cp vs x/c data could be stored in a data file FOROO9.DAT). 


This requires some knowledge of Fortran programming. 
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a. Sample Input to U2dif.f 
4 


ХЖХХХХХХАЖХХХХХХХЖХХХХХХХХХХХХХХХХХХХХХХХХХХХХХХЖХХЖХХХХЖХХЖХХЖХЖАХж жж 
AIRFOIL TYPE :NACAS0012 AIRFOILE 
NLOWER = 50 , NUPPER = 50 

ЖХХАХХХХХХХХХХХХХЖХХХХЖХХХХХЖХАХХХЖХХХХЖХХХХХХХЖХХХХХХХХХХХЖЖАХХЖХХХЖ 

J 5S0 50 

12 

OOO ESSA. 10.67.10. 0% 0.25, 0.00 0.00 

Оо, 0.00 

Jn ODIO, 0.0001, 0.00 

Шээгт222012115522224227 02 271 067 1007,да 283, 1554 

TTITLE 

INPUT TITLE (ITITLE - NO. LINE) 

IFLAG NLOWER NUPPER 

AIRFOIL TYPE 

ALPI DALP TCON FREQ PIVOT UGUST VGUST 

DELHX DELHY PHASE 

НЕ DIS О ТАРК 

ial ia2 1а3 1а4 1а5 1аб 1а7 1а8 1а9 1а10 
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b. Sample Output to U2diüf.f 


DATA FROM fort.l: (IRIS) OR FOROOINPAT K TERENE) 


* k k k k k kk KCk CkCKCK K koK KK CK KCk ACKCk k kCKCKCKCKCk KOXKUK KCKCAKCK K KCKCKCKCk KCKUKUKU KUX KCKCKUXCKCKCKCKCKUK UK KOK UK 


AIRFOIL TYPE : NACA OO12 AIRFOIL 
NLOWER = 50 , NUPPER = 50 


< k k k kk kk kk xx xxw xww W| xxx xww 


— — — —— — - — = — — — _  — — — — р ا ل‎ — — V oe ши 
— — — — е шш [l — —  — — D — —  — ` -—- -— — — — — — SS — — - | 





IFLAG (0:МАСА, l:INEUT) = 0 
NO. PANELS UPPER SURFACE = 50 
NO. PANELS LOWER SURFACE = 50 
INITIAL ANGLE OE ATTACK = 0.0000 
FINAL ANGLE OF ATTACK = 1525400 
RISE TIME шш O o 
REDUCED FREQ. FOR OSCIL = 0.00900 
PIVOT POINT = 0.2509 
STREAMWISE GUST VELOCITY = 0.0000 
PERPENDICULAR GUST VELOCITY = U 0000 
X AMPLITUDE OF TRANS OSCILL. = 0.0000 
Y AMPLITUDE OE TRANS 92 — 0709000 
PHASE OF TRANS OSCILL. = 0.0000 
FINAL TIME - 11.0000 
INITIAL TIME STEP FOR S.S. = ово 
TOLERANCE FOR CONVERGENCE = ПЛ! 
FACTOR BY WHICH DTS ADJUSTED = 0.0000 
COORDINATES OF AIRFOIL NODES 
Х/С Y/C 
1.000000 0.000000 
0.999013 -0.000141 
0.996057 -0.000562 
0.984292 0.002222 
0.991144 0.001258 
0.996057 0.000562 
0.999013 0.000141 
1.000000 0.000000 
AIRFOIL PERIMETER LENGTH = 2.039290 
STEADY FLOW SOLUTION AT ALPHA S 0.000000 
0], X(J) Od} GAMMA CP (J) V (J) 
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152433 
211236 
848815 
875229 
898071 
916416 
222360 
946548 
28293 
21179 


EU. 
e 
= ОР 
و‎ 
20, 
=e 
203 
x. 
ээ 
E 


22110 
"71173 
1259395 
‚946547 
89523959 
2216415 
2598070 
2596227 
. 848812 
511232 


дос 


. 752433 


.00000024 
. 00000025 
. 00000025 
.00000024 
. 00000024 


0.000001 


— 
— 


0.010000 


0143 
0036 


GAMMA 
-0:2209'67E-—06 
ОРАЗЕВОВЕ-05 
0.412307E-05 
(7211209Е-08 


1 0.999507 -0.112073 0.000000 0.433845 
2 0.997535 -0.120842 0.000000 0.341896 
3 0.993600 -0.126091 0.000000 0.279513 
4 0.987718 -0.129448 0.000000 0.232223 
5 0.979910 -0.131624 0.000000 0.193468 
6 0.970208 -0.132948 0.000000 0.160181 
7 0.958651 -0.133600 0.000000 0.130705 
8 0.945283 -0.133699 0.000000 0.104047 
9 0.930159 -0.133328 0.000000 0.079564 
10 0.913336 -0.132558 0.000000 0.056812 
90 0.894883 -0.131445 0.000000 0.035461 
91 0.913336 -0.132561 0.000000 0.056813 
92 0.930159 -0.133333 0.000000 0.079566 
93 0.945283 -0.133704 0.000000 0.104048 
94 0.958651 -0.133607 0.000000 0.130706 
95 0.970208 -0.132956 0.000000 0.160183 
96 0.979910 -0.131634 0.000000 0.193471 
97 0.987718 -0.129461 0.000000 0.232226 
98 0.993600 -0.126109 0.000000 0.279518 
99 0.997535 -0.120869 0.000000 0.341902 
100 0.999507 -0.112063 0.000000 0.433845 
COMPARISON OF G A MMA S 
GAMMA FROM KUTTA CONDITION: -0 
GAMMA FROM CONTOUR INTEGR (TRAIL EDGE): -0 
GAMMA FROM CONTOUR INTEGR (MIDPOINTS): -0 
GAMMA FROM BOX INTEGR (OFF THE CONTOUR): -0 
GAMMA FROM PRECISE CONTOUR ІМТЕС (6 РТ): -0 
CD = 0.000324 CL - -0.000000 cM = 
x k ХХХХХХХХХХХХХХХХХХЖЖХХХХХХХХХХ АХЖ 
ккк BEGIN UNSTEADY FLOW SOLUTION **** 
< < k xk x xxx xx ks x ж ж ss sk k kx sx kx 
TIME STEP TK = (02:02 0000 TK - TKMI 
ALPHA(T) - (73120021 ОМЕСА (Т) - -0.00 
U(T) - 0.000000 V (T) "7. 
NITR VXW VYW WAKE THETA 
0 1.000000 0.000000 0.010000 0.000000 
1 0.834730 0.000244 0.008347 0.000292 
2 0.827773 0.000255 0.008278 0.000309 
3 0.827452 0.000256 0.008275 0.000309 
CONVERGED SOLUTION OBTAINED AFTER NITR = 3 
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J X (J) Q(J) GAMMA EN) V (J) РЕЕОЖжВИТ (9) CURR PHI (@ 
1 0.999507 -0.107774 0.000004 2 222728 4 22. 0.032661 
2 0.997535 =0.117297 0.000004 02521094 2 117092 09125115 0.033045 
3 0.993600 -0.123077 0.000004 9.2952 02 849129: ЕЕ 0.055680 
4 0.987718 -0.126843 0.000004 1050292551 O 010011. 0.034404 
5 0.979910 -0.129345 0.000004 20 929409050 -0959990 921 221 0:095 2 Ыы 
б 0.970208 -0.130932 107000004 999061295 0790176627 QUIS UNS 0.056098 
7 0.958651 -0.131801 0.000004 9.12278 0593525 3] с 0703075 
8 0.945283 -0.132080 0.000004 05108051 059466DI ЕТ US 0.037463 
9 0.930159-0.131864” 02000004” 2129 672эЕ0 54220999 1-8 0.038041 
10 0.913336 -0.131225 -0:000002422202 029 5722:0 и 0.038482 
90 0.894883 -0.132662` 30. 09000860 222222:22г7 и 0.038794 
91 0.913336 -0.133893 0.000004 11025028 0227 7208 2 РО О Осо сма 
92 0.930159 -0:134798 02000004 КОКО 7710635 ЧО 5 o On ОБ 0.038099 
93 0.945283 -0.135322 0.000004 0.102044 0.946404 0.037474 0.037486 
94 0.958651 -0.135406 0.000004 9091029134 707232). 22 з 0.036724 
95 0.970208 -0.134972. 05000004 ОШ сг гое. 0.036025 
96 0.979910 -0.133914 000004 CA О 2202 Ег 0.03521 
97 0.987718 -0,132066 0:000004- 2022-2097 Ё и 0.034411 
98 0.993600 -0.129123 В ола 19879 ЗА ЕЕ EE 0.033664 
99 0.997535 =0.124414 02000004” 0 О 07 0 00 302 0.033044 
100 0.999507 -0.116362 9.000004 +0. 4 54960090: 75 UMEN OG 0.032639 
СОРТ (E GAMMAS 
GAMMA FROM KUTTA CONDITION: 0.00000411 
GAMMA FROM CONTOUR INTEGR (TRAIL EDGE): -0.00000165 
GAMMA FROM CONTOUR INTEGR (MIDPOINTS): 202090000139 
GAMMA FROM BOX INTEGR (OFF THE CONTOUR): 0.00000462 
GAMMA FROM PRECISE CONTOUR INTEG (6 PT): 07000003895 
CD = 0.000324 CL = 0. 0051056 CM = 090030519 
TRAILING VORTICES DATA 
M X (M) Y (M) CIRC 
17! 00413770 000018 0 000009 
TIME TER TK >. 0.020000 TK — THEE 0.010000 
ALPHA(T) - 0.000164 OMEGA(T) = -0.000285 
U(T) = 0.000000 V (T) = -0.000071 
NITR VXW VYW WAKE THETA GAMMA 
0 0.827452 0.000256 0.:008275 0 2 1-0 


170 


CONVERGED SOLUTION OBTAINED AFTER NITR - 0 


J X (J) Q(J) GAMMA CP (J) VT) PREV PHI (J) CURR PHI (J) 
ТШЕ тог 10102 2 00007171: 0.432928-0.753413 0.032661 0.032658 
27 


2 1:2-2-02212121729220:000011”-0.,341297  -0:8172045 09098045 0.033042 


.669285 -0.109493 0.000011 -0.121525 -1.061102 0.031441 0.031419 


20 0 

ОО 210212:272220702000тТ —0.138004 -1.069082 0.029361 22029339 
ПЕТ Т ТОП =0.155875 -1.077078 0.027013 02026992 
DEM v0 9099757] 0990000010: —0.1233227.—1.085108. 0.024398 0.024378 


| 


C. DIRECTIONS FOR THE EXECUTION OF INCOMPBL.F 
1. Steady State Case 
A. Execution of Incompbl.f is straight forward. First edit the input files, 


FOROO1.DAT and incompbl.dat, and set the desired flow conditions. 


TABLE 6. INCOMPBL.F INPUTS 


incompbl.dat IWAKE-Flag Wake Calulation 
NXT 
NW-No. of point in wake 
ITREND-No. of Iterations 
ITR(1)-Transition Flag Upper Surface 
ITR(2)-Transition Flag Upper Surface 
ISWPMX 
RL-Reynolds Number 
XCTR(1)-Transition Location Upper 
Surface 
IP-Print flag 


FOROOI.DAT Input AOA, pivot point and Airfoil 
Coordinates, and No. of Panels 


B.) Foraunix based machine the code is executed as follows: 

















incompbl «incompbl.dat» incompbl.out. 


TABLE 7. INCOMPBL.F OUPUT 


OUTPUT COMMENTS 





Note: The original Fortran code can be modified to store desired output in 
separate data files by inserting a WRITE (24,*) statement inside the required do-loop. 
(Example: Cp vs x/c data could be stored in a data file FORO24.DAT). This requires some 


knowledge of Fortran Programming. 


Ime 


a. Sample Input to Incompbl.f 
File Named fort.1 of FOROOI.DAT 
3 
С 
С NACA 0012 AIRFOIL 
С 
ALPI PIVOT 
4.000000 0.250000 
NLOWER NUPPER 
50 50 
х/С 
1.00000 0.98000 0.96000 0.94000 
0.88000 0.86000 0.84000 0.82000 
0.76000 0.74000 0.72000 22110000 
0,64000 062000 0.60000 9868000 
0552000 0.50000 0.48000 0.46000 
0.40000 0.38000 0.36000 0.34000 
0.28000 02220000 0,24000 8222000 
0.16000 0.14000 0.12000 ee 11000 
0.04000 0.02000 0.00000 0202000 
0.08000 0.10000 0.12000 0.14000 
0.20000 0822000 0.24000 0.26000 
9r 32000 0.34000 0.36000 0.38000 
0.44000 0.46000 0.48000 0.50000 
0.56000 0e 55000 0.60000 252000 
0.68000 me 70000 0 72000 0.74000 
0.80000 0.82000 S4000 0.86000 
0.92000 0.94000 0.96000 0.98000 
Y/C 
505126 20509205 0 006 74000938 
-0.01694 -0.01935 -0.02170 -0.02399 
-0.03056 -0.03264 -0.03467 -0.03664 
0.04222 P A 96T -0.04563 =0204723 
-0.05165 -0.05294 -0.05415 -0.05530 
-0.05803 -0.05868 -0.05923 -0.05966 
-0 05597 мэ 00753668 —0005911 8-025058 .38 
-0.05444 -0.05236 -0.04990 -0.04683 
—0.03231 2002382 0.00000 22177382 
0.04309 0.04683 0.04990 0.05236 
2 0252) 2305836 O. 05911 21052186 
0.05995 0.05966 0505923 0.05868 
0.05634 0205530 0, 05415 0.05294 
0.04878 0.04723 0.04563 0.04396 
0.03856 0.03664 0.03467 0.03264 
0.02623 0.02399 0702170 0.01935 
0.01196 0.00938 0.00674 0.00403 
File named incompbl.dat 
IWAKE NXT NW ITREND 
0 161 37 40 
ТТК (1) ЇТВ(2) ISWPMX RL 
0 0 1 3000000.0 
IP 
1 
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ROSES SO OE SOS (>) сә > 


zu 
=O: 
ЗГ 
er 
so 
=. 
=e 
=U 
7! 
‚05444 
2095997 
. 05803 
205165 
8222 
"5056 
.01694 
21901020249 


SOO Oe ae SOS 


„92000 
‚80000 
‚68000 
· 56000 
‚44000 
. 32000 
. 20000 
. 08000 
‚04000 
.16000 
. 28000 
. 40000 
292000 
. 64000 
. 76000 
. 88000 
. 00000 


ОС 
02623 
03856 
04878 
05634 
03995 
05737 
04309 


в!) 


0. 


10000 


CO SO GD Ss еее ее 


=O. 
227 
O 
E. 
t 
zu. 
U. 
=> 
.03842 
.05607 
.06006 
305120 
2550276 
.04042 
‚02842 
‚01448 


Sr СОО СО OS 


„90000 
23000 
. 66000 
. 54000 
.42000 
250000 
.18000 
.06000 
.06000 
.18000 
230000 
.42000 
‚54000 
. 66000 
. 18000 
.90000 


01448 
02842 
04042 
05026 
05726 
06006 
05607 
03842 


> 
e 
C 


1 


b. Sample Output from Incompbl.f 


NACA 0012 AIRFOIL 


INPUT DATA FOR INVISCID-FLOW CALCULATIONS 


ALPI- 4.0000 PIVOT- 072500 
NLOWER= 50 МОРРЕБ- 50 
COORDINATES OF THE BODY 
х/С 
1.000000 0.980000. 0.950000 3029400000909 20000 
0.880000 0.860000 0.840000 0.820000 0.800000 
0.760000 0.740000 0.720000 0.700000 0.680000 
0.640000 0.620000 0.600000 0.580000 0.560000 
0.520000 0.500000 0.480009 0.460000 0.440000 
0.400000 0.380000 0.360000 0.340000 0.320000 
0.280000 0.260000 0.240000 02220000 0110 
0.160000 0.140000 0.120000 0.100000 0.080000 
0.040000 0.020000 0.000000 0.020000 0.040000 
0.080000 0.100000 0.120000 0.140000 0.160000 
0.200000 0.220000 0.240000 0.260000 0.280000 
0.320000 0.340000 0.360000 0.380000 0.400000 
0.440000 0.460000 0.480000 0.500000 0.520000 
0.560000 0.580000 0.600000 0.620000 0.640000 
0.680000 0.700000 0.720000 0.740000 0.760000 
0.800000 0.820000 0.840000 0.860000 0.880000 
02920000 "0% 940000 0.960000 0.999090 T 7090090 
Y/C 
-0.001260 -02004030 5905006740 = ои 9580. 020 О 
-0.01694070019350 -030920200 17-0302 39 000076230 
-0.030560 -0:5032640029:0346 70 — 00365008. 000995 60 
-0.042220 -0.043960 -0.045630 -0.047230 -0.048780 
-0,051650 -0.052940 -0.054150 020553005 009555420 
20.058030 20.058680 -0.059230 207059660 ЕЕЕ 
30.059970 —058059660 -07099110 9059909» NNUS 7 0 
-0.054440 -0.052360 -0.049900 а) 046930, 090425090 
-0.032310 =0.023820 “0 000000” 09024220 ЕСЕТ 
0.043090 0.046830 0.049900 0.052360 0.054440 
0.057370 02058380 20205911072290596608 2097-2270 
0.059950 0:059660” 0 159230 08 0:36808 892120 
0,056340-20:055300:-10.054150т207:0529407Ш89 2 1550 
0.048780 0.047230 0.045630 0.043960 0.042220 
0.038560 0.036640 0.034670 0.032640 0.030560 
0.026230 0.023990 0. 021 70085 07295 0210 0 6941 
0.011960. 0.009380 0.00ё7:(07”9 0040 Яө 7 
ENVISCID FLOW RESULTS 
PANEL XP ҮР CE 
al 0.99000E*00 =O do 8 ЕР02 ) 255422500) 
2 0979005300 =0 448 T3E=02 Ота вео 
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CS OO OO See SS a е 


= Ok 
=; 
zu. 
=O. 
=; 
== О 
202 
=0 
1038420 
1056070 
. 060060 
2957260 
„050260 
‚040420 
. 028420 
. 014480 


соого Отоо 


. 900000 
. 780000 
. 660000 
. 240000 
. 420000 
. 300000 
.180000 
. 060000 
.060000 
.180000 
. 300000 
.420000 
. 540000 
. 660000 
. 780000 
.900000 


014480 
028420 
040420 
050260 
057260 
060060 
056070 
038420 


Ооо со ДЈ ду ла Од) 


ETE 
212 
ETS 
114 
115 
1: 6 
107 
PIS 
1816 


CL = 


OOO OO Oe aS 


CSS со отс лс СР) 


ооо C O O OQO O CO. COO CO С 


.95000Е+00 
. 93000Е+00 
.91000Е-00 
.89000Е-00 
au ODOESOO) 
2 001Е-00 
-52000Е-00 
2e POOCE TOO 


.79000Е-400 
.81000Е+00 
.83000Е+00 
.85000Е+00 
.87000Е+00 
.89000Е-00 
ОИ FOO 
.93000Е+00 
.95000Е-00 
.97000Е+00 
.99000Е+00 


ХР 


ВОО ЗЕ ЕО 
ОЕ О 
ООО 
Озак тот 
ОЕ От 
ПОРЕ НОЈ 
РЕЗЕ НОЛ 
“1 002Е-01 
:172:97Е-01 
.12543E+01 
ЗО! 
.14489Е+01 
159 ЈЕ ОЈ 
А РООРЕ 01 
:20759E-FOT 
посао 0: 
Не А о 
2598095501 
.38079E+01 


0.47937Е+00 


Oe OO © 2 о ого 


асосу о садад су ооо О СО 


. 4415E 02 0.1 ШТ ОЕ+0ОО 
IO 21E 01 0.9 AS E-Ol1 
NESUSGTE-OT 0. 74659E- 01 
Об БВ 06079 EO 
ILI EZ UL 0,49465Е-01 
.20502Е-08 0. 397995701 
222829Е-08 0:30916Е-01 
.25104Е-01 0. 218 Е-01 
21E UL -0.12424Е%00 
2506E 01 -0-105620Е 500 
но ASSET =0.88121E=01 
2ОЛО2Е-01 2212 2-1 
.18084Е-01 о 
. 1596642 0 =0.24583Е=01 
ПАРОВЕ OI 5995 05Е-02 
ОЕ 0:211598Е-01 
. 2966E OL 017984328501 
.43440Е-02 02159 146200 
.14480E-02 0S2 551L9E FOO 
INVISCID WAKE RESULTS 
YP ED 

wer 25-04 0.34783Е+00 
.93657Е-04 0.26481Е+00 
LOSSES 0.21845Е+00 
225-20Е-05 (22143 16Ет00 
SE OS OSS STOETOO 
T95 76CE=03 0221503-Е-00 
21 5101512)21448)7) 0:1208-17Е-00 
2232028502 25285550! 
Зоб Ел ОШ 26 ОБО 
2:1112Е-02 0 5740 65E-—01 
“50321502 0.44666Е-01 
"ШТО ШЕ ОРРВВАВЕ-0) 
ЩОБ Я БОЇ OSTEO 
2235 3ЇЕ-01 Рено ОВЕ ОТ 
219699 DE-—O 0 17 OTE OL 
во! U 20812502 
РОЗЕ о 2 114Е-02 
1059285009 ШЕП 5245-02 
ПЗА ТЭБ ВОО 0: 26592Е-—092 


INPUT DATA FOR BOUNDARY-LAYER CALCULATIONS 


IWAKE 


0 


T DEC) 


0 
TE 
1 


МХТ 
161 
ТТЕ?) 
0 


ДОМЕНА О“ "= 6 КБ 


NW 
217 


1 


175 


ITREND 


40 


3.00 


ACTR (1) 
0210 


o ooo o x 7K x x k CYCLE 


40 XC Ox x A x x 


SUMMARY OF THE DRAG, LIFT AND PITCHING MOMENT 
COEFFICIENTS WIIN THE COH 


CD,CL AND CM ARE EVALUATED FROM THE INTEGRATION OF CP 


CIC CE 


(ли ој N ٣۳ہ‎ 


(сол ее ie 


€ a M qeu cv aes E e 


CD 


2001993 
. 003184 
.002734 
. 002422 
2002221 


. 001684 
“O01751 
.001684 
200184 
‚001684 
200175 


€ OO Oo © 


SSO OO © 


CL 


‚001684 
‚ 001739 
. 001684 
. 004730 
. 001684 


‚427157 
‚427790 
.427165 
.427807 
.427167 
„427827 


(СО Со С 


СУ 5 


CM 


.427162 
.427756 
.427162 
"427770 
24271581 


.004733 
.004494 
.004731 
.004491 
‚004731 
‚004488 
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